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Abstract  / 

The  blanket  objective  of, the  four  tasks  of  the  research 

/ 

proposed  and  completed  under' this  contract  was  to  tailor  TRAPS 
(TRansient  Acoustic  Processing  System)  to  the  seismic  event 
monitoring  application  thus  providing  a  powerful  tool  for  the 
nuclear  monitoring  problem.  TRAPS  is  a  three  stage  process. 

Stage  One  is  mathematically  designed  to  perform  noise  elimination 
on  each  receiver  component.  -The  output, of  Stage  One  is  then 
beamformed  and  input  into  Stage  Two,  which  performs  signal 
enhancement.  -The'*  information  resulting  from  Stage  One/Two 
is  used  in  Stage  Three  of  TRAPS'  to  calculate  the  location  and 
classification  of  the  detected  event..  Therefore,  the  achieved 
goals  of  this  contract's  research  will  be  detailed  in  this  report 
by  relating  the  refinements  to  these  stages,  resulting  from  the 
contract  analysis. 

While  the  original  TRAPS,  as  pictured  in  Figure  1>  succeeded 
in  doubling  the  range  of  detection^  in  a  previous  seismic  study 
( 1 ) ,  it  did  so  with  a  standard  beamformer.  Task  1  was  th us> to 
ascertain  if  an  adaptive  beamformer  (AB)  incorporated  into  TRAPS 
would  produce  any  gains.  Task  1  research  efforts- determined 
TRAPS'  Stage  One  followed  by  a  standard  beamformer  was  superior 
to  an  AB,  yet  an  AB  strategically  placed  between  TRAPS'  Stage  One 
and  Stage  Two  did  provide} some  processing  gain. 

The  goal  of  Task  2  was  to  investigate  and  incorporate  into 
TRAPS  the  relevant  background  noise  statistics  for  increased 
TRAPS'  efficacy.  Task  2  research  efforts  did  result  in  the 
formulation  of  algorithm  refinements  to  TRAPS.  One  such  refine¬ 
ment  to  TRAPS  was  the  upgrade  of  the  original  adaptive  auto¬ 
regressive  (AR)  model  of  TRAPS  to  a  unique  cascaded  adaptive 
ARMA  (CAARMA)  model,  which  provided  significant  improvement  to 


TRAPS.  Anti  jamming  strategies  to  prevent  deterioration  to 
TRAPS  by  local  interference,  etc.,  were  then  included  into  TRAPS 
by  way  of  a  Slaved  Stage  One.  Additionally,  the  order  and 
adaptation  time  of  Stage  One  and  Two  of  TRAPS  was  varied  to 
improve  TRAPS'  ability  to  remove  bias  in  the  data. 

The  Task  3  goal  was  to  decrease  the  false  alarm  rate  of 
TRAPS  and  to  enhance  FK  space  definition.  It  was  proposed  to 
accomplish  Task  3  goals  by  means  of  a  signature  analysis  (SA) . 
Since  the  coefficients,  poles  and  zeros,  etc.  generated  by  the 
CAARMA  process  completely  characterize  the  shape  of  the  power 
spectrum  and  thus  an  event,  this  SA  characterization  mode  was 
composed  of  the  weighted  prediction  error  residual  (WPR)  of  each 
event  in  combination  with  the  event's  weighted  prediction  co¬ 
efficients  (WPC) .  The  SA  allowed  the  achievement  of  several 
advancements  to  TRAPS.  First,  full  waveform  or  total  energy 
was  appropriately  classified  into  its  various  wave  components. 
This  allowed  utilization  of  all  available  information  of  an  event 
for  improved  detection.  Once  an  event  was  identified,  the  poles 
of  the  CAARMA  model  were  uniquely  applied  to  estimation  of  time 
delay  among  phones.  The  resulting  estimation  was  input  to  a 
feedback  maximum  likelihood  scheme.  The  culmination  of  these 
research  efforts  produced  a  powerful  processing  system,  TRAPS  II. 

Finally,  Task  4  examined  TRAPS  for  its  possibility  of  a 
compact  digital  processing  system.  Analysis  indicated  TRAPS 
and  TRAPS  I  (see  Figures  1  and  2)  were  capable  of  real  time 
processing.  TRAPS  II,  as  depicted  in  Figure  3,  would  require 
an  array  processor  for  real  time  capability. 


1 . 0  SUMMARY 


TRAPS  (TRansient  Acoustic  Processing  System)  is  a  three 
stage  process.  Stage  One  is  a  hybrid  adaptive  filter,  which 
is  mathematically  designed  to  perform  noise  elimination  on  each 
phone.  The  output  of  Stage  One  is  then  beamformed  and  input 
into  Stage  Two,  another  hybrid  adaptive  filter  which  serves  to 
enhance  the  signal.  The  information  resulting  from  Stage  One/ 
Two  is  used  in  Stage  Three  of  TRAPS  to  calculate  the  range, 
bearing  and  classification  of  the  detected  event.  The  blanket 
objective  of  the  research  proposed  and  completed  under  this 
contract  was  to  improve  the  existing  detection  and  localization 
capability  of  TRAPS  by  means  of  tailoring  each  of  TRAPS'  three 
Stages  to  the  seismic  event  monitoring  application.  This  broad 
objective  was  broken  down  into  4  specific  Tasks. 

While  the  original  TRAPS  succeeded  in  doubling  the  range 
of  detection  in  a  previous  seismic  study  (1) ,  it  did  so  with  a 
standard  beamformer,  as  opposed  to  an  adaptive  beamformer  (AB) . 
Task  1  was  thus  to  ascertain  the  most  effective  mode  of  AB 
incorporation  into  TRAPS  and  specify  any  gains  resulting  from 
the  incorporation.  Task  1  research  efforts  determined  TRAPS' 
Stage  One  followed  by  a  standard  beamformer  was  superior  to  an 
AB,  yet  an  AB  strategically  placed  between  TRAPS'  Stage  One  and 
Stage  Two  did  provide  some  processing  gain.  Section  2  compares 
the  effectiveness  of  the  original  TRAPS  to  TRAPS  with  an  AB  and 
to  other  spectral  estimators  as  illustrated  by  means  of  a 


normalized  frequency  procedure. 

The  goal  of  Task  2  was  to  analyze  the  environment  within 
which  TRAPS  was  to  function.  The  knowledge  of  the  relevant 
background  noise  statistics  accumulated  by  this  study  would 
then  be  incorporated  into  TRAPS  for  its  improvement.  Task  2 
efforts  did  inspire  the  formulation  of  two  algorithm  refine¬ 
ments  to  TRAPS.  The  first  refinement  addressed  the  signal  can¬ 
cellation  phenomena  which  is  a  consequence  of  adaptive  inter¬ 
action  between  the  desired  signal  and  interference,  certain 
possible  false  signals  generated  that  defeat  or  partially 
deteriorate  the  adaptive  filter,  etc.  As  a  result,  a  sophis¬ 
ticated  slaved  algorithm  was  built  into  Stage  One  of  TRAPS  as 
an  anti-jammer  to  prevent  local  interference,  etc.  The  second 
change  to  TRAPS  was  influenced  by  the  knowledge  that  the  auto¬ 
regressive  moving  average  (ARMA)  model  often  provides  better 
spectral  estimates  than  the  more  specialized  AR  models.  Thus 
the  original  adaptive  AR  model  of  TRAPS  was  transformed  to  an 
adaptive  ARMA  model,  which  could  provide  source  parameterization 
improvement  in  a  low  SNR  environment,  while  maintaining  real 
time  processing  criteria.  These  refinements  are  presented  in 
Section  3.2.  The  resulting  model  was  then  upgraded  by  a  unique 
cascade  approach  to  the  problem.  This  cascaded  adaptive  ARMA 
(CAARMA)  model,  which  is  detailed  in  Section  3.3  and  3.4  pro¬ 
vided  excellent  results.  Additionally,  the  order  and  adaptation 
time  of  Stage  One  and  Two  of  TRAPS  was  varied  to  improve  its 
ability  to  remove  various  bias  in  the  data,  as  discussed  in 
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Section  3.5.  The  end  result  was  a  processing  system  that  had 
excellent  results,  as  presented  in  Section  3.6. 

The  goal  of  Task  3  was  to  decrease  the  false  alarm  rate 
of  TRAPS  and  to  enhance  FK  space  definition  by  means  of  a 
signal  characterization  mode.  In  order  to  demonstrate  the 
logic  for  the  direction  taken  in  Task  3  it  should  be  noted  that 
the  coefficients,  poles  and  zeros,  etc.,  generated  by  the  CAARMA 
process  completely  characterize  the  shape  of  the  power  spectrum 
and  thus  an  event.  Thus  the  similarity  of  one  event  to  another 
can  be  verified  by  comparing  the  parameters  which  result  from 
Stage  One  and  Two  of  TRAPS.  The  signature  of  one  event  was 
then  compared  to  the  signature  of  another  event  by  means  of  the 
weighted  prediction  error  residual  (WPR)  of  each  event.  These 
WPR's  in  combination  with  the  weighted  prediction  coefficients 
(WPC)  served  as  the  foundation  for  a  statistical/structural 
classification  system,  or  pattern  recognition  algorithm  (PRA), 
as  presented  in  Section  4.3. 

The  PRA  allowed  the  achievement  of  several  improvements  to 
TRAPS.  First,  such  a  PRA  is  able  to  discriminate  between 
various  types  of  explosions  and  earthquakes.  Also,  the  infor¬ 
mation  resulting  from  the  PRA  isolated  and  then  classified  the 
various  waves  of  an  event.  This  PRA  was  also  capable  of  identi¬ 
fying  phase  irregularities  that  characterized  specific  field 
data.  This  technique  uses  the  full  waveform  or  total  energy, 
thus  allowing  utilization  of  all  available  information  of  an 
event  for  determination  of  detection,  location  and  classification. 


There  is  no  penalty  built  into  the  system  for  "missing"  infor¬ 
mation.  The  advantage  lies  in  the  fact  that  whatever  informa¬ 
tion  is  present  is  appropriately  weighted  to  provide  maximum 
efficacy.  Once  an  event  was  identified  the  poles  of  the 
CAARMA  model  were  uniquely  used  for  time  delay  estimation. 

This  innovative  approach  was  initiated  by  earlier  results  that 
the  nature  of  the  pole  migration  served  as  a  robust  event  de¬ 
tector.  This  time  delay  estimation  was  input  to  a  feedback 
maximum  likelihood  scheme.  The  culmination  of  these  research 
efforts,  which  produced  a  powerful  processing  system,  is  dis¬ 
cussed  in  Section  4.4. 

Finally,  Task  4  examined  TRAPS  for  incorporation  into  a 
compact  digital  system.  Once  the  reference  sets  of  the  PRA  were 
established,  it  was  found  that  human  intervention  was  not 
necessary.  The  only  limit  to  the  TRAPS  resulting  from  the 
above  research  was  that  real  time  processing  could  occur  only 
with  the  aid  of  an  array  processor. 
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2.0  TASK  1:  EVALUATION  OF  TRANSIENT  ACOUSTIC  PROCESSING 


SYSTEM  (TRAPS')  FOR  PROCESSING  GAINS  DUE  TO  AN  ADDED 
ADAPTIVE  BEAMFORMER 


2 . 1  Introduction  to  the  Fundamental  Hybrid  Model  of 

TRAPS 

Section  2.1  shall  present  the  mathematical  founda¬ 
tion  of  the  original  TRAPS  and  its  relationship  to  the  physics 
of  the  nuclear  monitoring  problem.  This  process  will  enable 
the  reader  to  then  comprehend  the  reason  for  TRAPS*  initial 
success  in  seismic  event  detection  and  location  (1) ,  the  logic 
behind  the  various  improvements  to  TRAPS  and  finally  insight 
into  the  unique  capability  and  consequential  success  of  the 
improved  TRAPS,  TRAPS  I  for  the  nuclear  monitoring  application. 
It  is  important  to  understand  that  TRAPS  is  in  essence  a  three 
stage  algorithm.  These  stages  are  (as  illustrated  in  Figure 
1): 

*  Stage  One  -  Noise  Elimination 

*  Stage  Two  -  Signal  Enhancement 

*  Stage  Three  -  Range  and  Bearing  Estimation 

This  report  shall  first  concentrate  on  some  of  the 
terminology  and  theory  of  the  adaptive  autoregressive  (AR) 
model  which  Stage  One  and  Two  of  the  original  TRAPS  was  based. 
Let  x(n)  «*  x(tn)  °  x(nAt);  n  =  0,  1,  2,...,  denote  a  particular 
seismic  wave  component  at  the  geophone  location.  x(n)  is  then 


considered  to  have  originated  from  a  source  located  some  un¬ 
known  distance  from  the  geophones.  The  received  sequence  is 
assumed  to  satisfy  a  difference  e^ nation  of  the  sort: 

°  0 i ^ t - 1  +  ^  2  ^  t  -  2  +  *’*  +  4p  ^  t - p  +^t  ^ .  1 

where  the  0  are  constants,  and  the  forcing  function  E  is  a 
P  ^ 

stochastic  variable  with  Gaussian  statistics  during  an  interval 

T  i  nat  i  T.  and  is  zero  for  nAt  outside  this  interval.  Given 
o  i 

that  (To,  Tj)  =  (-—  ,-)  and  that  £t  is  stationary  this  model 
would  be  called  an  autoregressive  process  of  order  p  (2) . 

For  this  adaptive  AR  model  of  TRAPS  the  finite  dura¬ 
tion  stochastic  function,  £t,  can  be  considered  to  represent 
the  unstructured  energy  in  the  wave,  and  the  coefficients, 

0p ,  can  be  thought  of  as  representing  the  filtering  effect  of 
the  earth  and  the  initial  source  structure.  The  received  wave 
is  modeled  as  the  output  of  a  linear  recursive  system  which  at 
some  time  Tq  is  excited  by  Gaussian  noise,  and  at  some  time  T^ 
is  de-excited.  While  this  is  still  an  autoregressive  process 
of  order  p  (3),  the  sequences  of  the  same  type  as  eqn  (2.1)  in 
hybrid  adaptive  form  is  referred  to  as  "transient  autoregressive. 
The  order  to  the  transient  autoregressive  process  is  critical 
and  various  methods  of  selection  have  been  studied  extensively 
(4-7).  The  method  chosen  for  identifying  the  order  of  the  pro¬ 
cess  and  the  order  we  chose  for  nuclear  monitoring  is  addressed 
in  Section  3.5.  This  allows  us  to  consider  p  as  known  for  this 


V 


discussion.  However,  we  do  not  know  either  [ 0p]  or  [£tl  and 
in  order  to  circumvent  this  difficulty  we  employ  a  hybrid 
adaptive  filtration  technique.  We  summarize  here  the  adaptive 
least  mean  square  (LMS)  filter  originally  used  in  these  data 
studies.  The  theoretical  basis  of  the  adaptive  LMS  algorithm 
is  an  extension  of  work  by  Widrow  (8)  and  Griffiths  (9),  which 
can  be  referenced  for  further  detail. 

The  adaptive  filter  constructs  the  sequence  of  pre¬ 
diction  filter  coefficients  [0p(n)]  which  converge  to  true  co¬ 
efficients  as  n-»«°,  by  first  predicting  a  input  2(n)  at  time  n 
from  a  linear  combination  of  previous  values  of  the  actual 
input  Z(n).  Thus, 


ZOO  -  E  0_(n)Z(n-p) 
p=l  p 


This  allows  an  error  sequence  E(n) , 


2.2 


6  00  =  Z(n)  -  Z(n) 


2.3 


to  be  calculated.  The  set  of  filter  coefficients  0-.*,0  *,...,0  * 

X  z  p 

2 

which  yields  the  minimum  mean-square-error  E[&  (n)]  is 


E[e2(n)]  =  E[[Z(n)-£(n)]2]. 


2.4 


This  equation  is  a  quadratic  function  of  the  filter  coeffi¬ 
cients  and  has  a  unique  minimum  value  which  is  achieved  by 


t 
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[ 0p* ] CIO) •  The  coefficients  are  then  given  by  the  LMS 
algorithm 


0p(n+l)  -  0p(n)  +yu[Z(n) -Z(n) ]  Z(n-p)  2.5 

where  /x  is  the  proportionality  constant.  The  choice  of  yu. 
determines  both  the  rate  of  convergence  of  the  processor  and 
the  steadystate  mean-square-error.  Convergence  of  the  algo¬ 
rithm  is  assured  for  stationary  autoregressive  data  [Z(n)]  of 
order  P  provided 

I 

- 

ju.  *=  P  2.6 

whereof  is  the  input  power  level  and  0<4<2.  The  adaptation 
time  constant  of  the  algorithm  expressed  in  data  points  is 

t  -  — — — r  •  2-7 

Ln(ly.rzn 

The  method  of  selecting  the  adaptation  time  as  well  as  its 
effect  on  the  processing  gain  will  be  presented  in  Sections 
2 . 2  and  3.5. 

The  1’irst  Stage  of  TRAPS  uses  this  adaptive  tran¬ 
sient  AR  model  to  filter  out  the  ambient  noise.  The  ambient 
noise  field  is  considered  to  be  the  sum  of  a  number  of 
travelling  waves.  The  sampled  time  trace  at  a  fixed  spatial 
point  of  each  of  these  waves  is  assumed  to  be  a  transient 
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autoregressive  sequence  as  defined  earlier.  Therefore,  the 
measured  time  trace  of  the  ambient  noise  field  at  a  fixed 
point  is  a  sequence  which  is  the  sum  of  a  number  of  transient 
autoregressive  sequences.  Such  a  sequence  is  not  in  general 
autoregressive.  The  following  assumptions  are  made  regarding 
the  ambient  noise  field: 

1)  there  is  a  dominant  (energy  wise)  AR  component 

2)  this  component  is  stationary  or  at  worst  slowly 
varying 

3)  an  adaptive  least  mean  square  (LMS)  filter  with 
appropriate  parameters  converges  to  (and  tracks) 
this  component. 

These  assumptions  are  justified  empirically  by  a 
large  number  of  ambient  spectra  that  we  derived  from  the  NORESS 
data.  The  assumptions  are,  moreover,  justified  physically  by 
virtue  of  the  bandpass  property,  with  respect  to  wavenumber, 
of  a  layered  half -space.  The  bandpass  characteristic  of  the 
earth  channel  leads  to  dominant  energy  frequencies  as  well 
as  a  temporal  extension  of  signals.  The  superposition  of  a 
number  of  such  signals  results  in  a  stationary  signal  as  in 
the  classical  analysis  of  filtered  shot  noise  (11). 

The  original  TRAPS  algorithm  presented  in  this 
section  was  designed  to  model  all  distinct  signals  which  are 
measured  by  the  geophone  array.  All  wave  activity  travelling 
the  earth  originates  from  some  source,  and  therefore  in  this 
sense  all  measurements  are  due  to  sources  of  the  type  we  are 


interested  in.  However,  as  each  signal's  distance  from  its 
origin  increases  differential  absorption,  normal  spreading 
and  viscous  effects  cause  signals  to  overlap  one  another  more 
and  are  generally  of  lower  individual  amplitude  and  jointly 
of  less  pronounced  structure.  When  these  signals  are  received, 
the  net  intensity  may  be  significant.  Stage  One  of  TRAPS 
attempts  to  eliminate  as  much  of  this  quasi-stationary  back¬ 
ground  as  possible  by  means  of  a  particular  adaptive  filter. 

Let  us  suppose  that  we  have  knowledge  of  the  joint 
ambient  autoregressive  equation  postulated  in  item  1  above. 

We  remark  that  this  equation  can  be  obtained  adaptively  as 
proposed  in  item  3  above  or  by  means  of  an  asymptotic  surface 
wave  analysis  based  upon  a  detailed  knowledge  of  the  geographi¬ 
cal  layering  at  the  test  site  (12) .  The  measurement  may  then 
be  written  as 

x(n)  =  z(n)  +  u(n)  2.8 

where  u(n)  is  the  dominant  ambient  component  so  that  u(n) 
satisfies  the  "ambient  equation": 

L  £u(n)]  =  w(n)  2.9 

and  z(n)  is  the  residual  component  not  satisfying  this  equation 
(13).  The  notation  L[u(n)j  denotes  the  autoregressive  operator 
L  acting  on  the  vector  of  data  values  u(n)  known  at  time  n. 


The  variable  w  is  the  excitation  variable,  which  by  assumption 
2  above  is  essentially  white  noise. 


The  noise  cancellation  proceeds  as  follows: 


Define 


£(n)  -  L  [x  (n)] 


2.10 


Then  we  see 


£ (n)  -  y (n)  +  w(n) 


2.11 


where  y (n)=L  [ Z (n)]  is  transient  autoregressive.  That  y(n) 
is  transient  autoregressive  follows  from  the  commutativity  of 
constant  coefficient  difference  operators.  The  earlier  defini¬ 
tion  of  [Tq,  T^ }  shall  be  extended  very  slightly  to  allow  the 
interval  [Tq,  Tjl  to  be  a  union  of  such  intervals.  The  defini¬ 
tion  itself  does  not  require  the  noise  to  be  homostochastic. 

For  these  LMS  filters  acting  on  the  individual 
phones  the  adaptation  parameters  must  be  chosen  such  that 


TA&rat»  Tg  2.12 


where  TA  is  the  ambient  time  scale  and  Tg  is  the  signal 
duration.  This  time  scale  allows  TRAPS  to  track  and  then 
eliminate  the  ambient  noise  structure  without  disturbing  the 
transient  signal  structure.  In  summary,  Stage  One  gives  us  a 
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remarkable  (and  computationally  fast)  SNR  improvement,  due  to 
the  following  virtues: 


1)  Independent  adaptive  filters  implemented  for 
each  geophone,  K,  (and  channel)  in  an  array 
will  converge  to  the  ambient  operator  which 
is  local  to  that  geophone.  The  resulting  noise 
sequences  w^fn)  will  be  independent. 

2)  The  sequences  E^fn)  can  then  be  combined  in  a 
manner  chosen  to  accentuate  the  spatially 
coherent  structure  of  any  incoming  signals 
which  are  rapidly  varying  with  respect  to  the 
ambient  time  scale.  The  independence  of  the 
sequences  w^fn)  will  produce  maximum  signal  to 
noise  ratio. 

Once  the  first  Stage  of  TRAPS  has  performed  the 
necessary  noise  elimination,  the  filtered  data  is  then  pro¬ 
cessed  by  Stage  Two  of  TRAPS.  Stage  Two  gives  mathemetical 
priority  to  signal  enhancement.  In  order  to  describe  an  in¬ 
coming  wave  of  interest  we  let  yp(n)  denote  the  observed  geo¬ 
phone  measurement  on  geophone  p.  Then 


z(n)  s  E  yn(n)  cos  (k.x  +  k?y  ) 
p=l  P  1  P  *  P 


2.13 


where  (Xp,y^)  are  the  coordinates  of  the  p**1  geophone  and 


kj  ■  _2tf  cos  0 


2.14 
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2TT  sin  0 

X 
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where  X  and  0  are  wavelength  and  bearing  of  an  incoming  plane 
wave,  respectively.  The  instantaneous  power  spectrum  of  the 
detected  wave, 

H  (f;n)  of  z(n)  at  time  t  ■  nat  and  frequency  f, 
z  n 

is  given  by 


^£;n) 


[h(f;n)]2 


M 

h(f;n)  =  1  -  am(n)  exp  (-i2tifmat) 


2.16 


2.17 


where  am(n)  are  adaptive  coefficients  modeling  z(n)  as 

„  M 

Z (n)  =  £  a  (n)  Z(n-m)  2.18 

m=l  m 

am(n+l)  ■  am^n^  +  UZ (n-m) [z (n)  -  z(n)]  2.19 


For  the  filters  acting  on  this  beamformed  data  the 
adaptation  time,  n,  varied  from  one  event  to  another  in  order 
to  incorporate  the  local  variance  of  each  event. 

The  presence  of  a  signal  was  determined  by  sharp 
changes  in  the  value  of  the  spectrum  (2.16)  within  a  frequency 
range  acceptable  to  physical  laws  that  govern  the  nuclear 
monitoring  situation.  Filtration  in  this  narrow  band 
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successfully  detected  the  provided  events.  The  range  of 
detection  for  previous  seismic  applications  (1)  was  more  than 
doubled. 

Figures  4  and  5  are  provided  to  illustrate  just  how 
effective  is  the  first  stage  and  second  stage  of  TRAPS  for  the 
NORESS  application.  Figure  4  is  the  output  of  just  Stage  Two 
of  TRAPS.  Figure  5  is  the  output  achieved  by  using  Stage  One 
and  Stage  Two  of  TRAPS. 


2.2  Mathematical  Detail  of  Task  1  Evaluation  Procedure 


Task  1  of  the  Contract  was  to  expand  TRAPS  to  in¬ 
clude  an  adaptive  beamformer  (AB) .  It  was  necessary  to 
evaluate  the  positioning  of  the  adaptive  beamformer  into  TRAPS 
so  as  to  determine  the  position  which  provided  greatest  pro¬ 
cessing  gain.  The  performance  criterion  used  to  judge  the 
processing  gain  of  the  various  positions  an  AB  in  TRAPS  was 
a  two-sinusoid  frequency  resolution  (14).  This  method  may  be 
best  understood  by  considering  two  sinusoids  with  frequencies, 
f1  and  f2.  The  center  frequency  of  these  two  frequencies, 
f  ,  is  (f  +  f_)/2.  The  condition  for  two-sinusoid  resolution 

C  1  fa 

in  this  evaluation  is  defined  as  the  frequency  separation, 

&f  *  [f1  -  f2]  at  which  the  power  spectral  density  evaluated 
at  the  center  frequency,  S(fc),  is  equal  to  the  average  of  the 
power  spectral  densities  evaluated  at  the  two  sinusoid  fre¬ 
quencies  ,  i .  e .  , 


S(fc)  =  2[s(f1)  +  S(f2)] 


2.20 


This  condition  is  demonstrated  in  Figure  6.c.  The  motivation 
for  the  definition  of  the  condition  for  resolution  was  to  pro¬ 
vide  a  common  method  that  could  be  applied  to  the  various 
spectral  estimation  under  investigation.  The  definition  is 
independent  of  any  property  of  a  specific  spectral  estimation 
method.  It  is  also  computationally  simple  to  evaluate.  Eq. 


2.2Q  does  not  require  exact  knowledge  of  the  location  of  the 
spectral  peak.  Also,  Eq.  2.20  was  beneficial  in  that  the 
just-resolved  condition  of  Figure  6.c.  occurs  at  a  point  near 


the  resolution  limit,  so  that  a  slight  shift  of  f^  and  f2 
toward  f  would  merge  the  peaks  into  one. 

The  measure  of  frequency  resolution  is  thus  provided 
by  a  dimensionless  quantity  R,  or  normalized  resolution, 

R  =  2tr  a  taf  2.21 

where  At  is  the  sampling  interval  in  seconds  and  A f  is  the 
frequency  separation  in  Hz  point  resolved  according  to  Eq.  2.20. 
This  was  viewed  to  be  an  optimum  definition,  since  the  only 
frequency  variable  required  in  the  computer  simulations  for  all 
spectral  estimation  techniques  was  the  radian  frequency  X, 
where 

9 

^  -  2trfAt  -  X<1r  2.22 

and  f  is  the  true  frequency  in  Hz.  Thus,  it  was  unnecessary 
to  specify  an  actual  sample  interval  At  or  frequency  f  in 
order  to  run  the  simulations.  The  condition  of  resolution 
occurs  when 


S 


( 


li  *  h 


2 


2 


[S(\)  +  stt2)] 


2.23 


The  true  resolution  in  Hz  is  obtained  from  the  normalized 
resolution  by 

Af  -  R/2trAt  2.24 


For  the  following  comparisons  use  was  made  of  a 
modified  normalized  resolution 

Rn  =  2trNAtAf  2.25 

plotted  as  a  function  of  N,  the  number  of  data  points,  or  auto¬ 
correlation  lags.  Given  the  sample  interval.  At  and  N,  plus 
the  value  of  RN  determined  for  the  spectral  estimation  technique 
being  evaluated,  the  frequency  resolution  in  Hz  will  be  given 
by 


Af  =  R../2frNAt  2.26 

N 

Note  that  smaller  values  of  RN  imply  better  frequency  resolution. 

The  effect  of  noise  on  the  resolution  of  the  various 
spectral  methods  evaluated  was  measured  by  running  an  ensemble 
of  computer  simulations  that  added  white  noise  to  the  two 
sinusoids.  During  each  resolution  measurement,  the  noise  was 
held  constant  while  the  frequency  separation  of  the  sinusoids 
was  gradually  decreased  until  the  resolution  criterion  was  met. 
For  each  series  of  measurements,  this  procedure  was  repeated 
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over  independent  noise  epochs.  The  initial  phases  of  both 
sinusoids  were  fixed  at  0°.  The  center  frequency  between  the 
sinusoids  was  at  half  the  foldover  frequency,  or  the  center 
of  the  Nyquist  bandwidth. 

The  goal  of  Task  1,  as  mentioned  earlier,  was  to 
evaluate  whether  or  not  there  was  a  differential  gain  of 
TRAPS  once  expanded  to  include  an  adaptive  beamformer  (AB) . 
Evaluation  needed  to  also  demonstrate  the  various  positions  of 
an  AB  within  TRAPS  (front  end,  between  Stage  One  and  Two  of 
TRAPS,  etc.)  and  their  respective  differential  gains.  Adaptive 
beamforming  in  the  time  domain  can  basically  be  categorized 
into  two  mathematical  concepts,  the  constrained  minimum  power 
and  the  unconstrained  minimum  "erroT."  Due  to  this  fact  two 
versions  of  TRAPS  were  prepared  for  AB  implantation.  The  first 
algorithm  was  based  on  Wen-Wu  Shen ' s  "A  constrained  minimum 
power  adaptive  beamformer  with  time-varying  adaption  rate"  (15) , 

9 

which  we  shall  hereafter  refer  to  as  Shen's  algorithm.  Shen's 
algorithm  was  chosen  due  to  its  ease  of  implementation  into 
TRAPS  and  its  reported  success  with  seismic  data.  The  second 
algorithm,  which  fit  into  the  unconstrained  minimum  error 
category,  was  based  on  an  algorithm  recently  developed  by 
Bernard  Widrow  and  Kenneth  Duvall  at  Stanford  University  (16) . 
This  algorithm,  hereafter  referred  to  as  Widrow’ s  algorithm, 
was  implemented  into  TRAPS  for  comparison  also  due  to  its  re¬ 
ported  success  in  seismic  applications.  The  reader  is  referred 
to  Shen’s  (15)  and  Widrow's  (16,  17)  articles  for  mathematical 
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detail  of  each  AB,  since  the  intent  of  this  report  is  not  to 
duplicate  well  written  previously  published  algorithms. 

To  make  the  Task  1  evaluation,  it  was  believed  that 
a  standard  power  spectral  density  method  should  be  included  in 
the  analysis  to  serve  as  a  benchmark.  A  form  of  the  periodo- 
gram  was  chosen  since  it  is  a  standardized  method  of  power 
spectral  density  (18)  .  The  modern  form  of  the  periodogram  is 
computed  as  the  modulus  of  a  modified  Fourier  transform. 

a  ^  N  - 1 

Sper(f)  =  rj—  [*t  £  Xn  exp(-j2frfnAt)  ]2  2.27 

n=0  n 

for 

n  <  f  <  f 

'  1  '  rNyquist 

where 

^Nyquist  =  t*ie  folding,  or  Nyquist  frequency 

At  =  sampling  interval 
N  =  number  of  data  samples 
NAt  =  total  observation  time 

=  x(nAt),  the  data  samples  for  n=0,...,  N-l. 

The  function  described  by  Eq.  2.27  is  a  cyclic  function  of  f, 
repeating  its  values  at  intervals  of  1/AtHz. 

By  evaluating  Eq.  2.27  at  only  a  set  of  frequencies 
spaced  equally  apart  by  1/NAt,  the  DFT  method  of  spectrum 


♦ 
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analysis  is  obtained: 


N-l 

F  -  At  T~  x  exp  (-j2lrnm/N)  2.28 

ro  n-0  n 

§nFT(m)  -  1  IF  I2,  m  ■  0 , . . . ,  H-l  2.29 

DFT  RZt  m 

The  exponential  argument  is  not  a  function  of  the  sample 
interval.  At,  due  to  the  particular  selection  of  the  frequency 
spacing.  That  is,  the  same  number  of  transform  "samples"  as 
data  samples  are  computed.  If  N  is  chosen  to  be  2  ,  for  r  some 
positive  integer,  then  a  computationally  efficient  algorithm, 
known  as  the  FFT,  may  be  used  to  solve  for  the  Fm  of  Eq.  2.28. 


2.3  Technical  Results 


Figure  7  is  the  plot  of  the  normalized  resolution 

R  versus  M  for  mean  resolution  of  a  standard  FFT  (periodogram) 

M 

and  the  various  TRAPS  algorithms  resolution  at  three  different 

SNRs .  The  SNR  is  defined  as  a  function  of  the  sinusoid  power 

2 

Ps  and  the  noise  variance,  =CTn  » 

SNR (in  dB)  =  10  log  Pg/PN  2.30 

A  fact  to  be  noted  from  Figure  7  is  that  the  mean 
resolution  of  the  periodogram  is  independent  of  the  SNR.  How¬ 
ever,  the  variance  about  the  mean  increases  with  decreasing 
SNR.  Table  1  demonstrates  this  phenomena  for  N  ■  256.  The 
SNR  independence  is  not  surprising  since  the  windowing  is  the 
key  contributing  factor  that  limits  the  periodogram  PSD 
estimate  resolution. 

The  AR  PSD  resolution  is  also  a  function  of  the  SNR. 
It  can  be  seen  from  Figure  7  that  the  AR  technique  produces 
significantly  higher  resolution  than  the  conventional  technique. 
For  example,  in  Figure  7  for  N  =  33  and  20db  SNR,  there  is  a 
significant  improvement  in  resolution  using  the  AR  technique 
over  the  averaged  periodogram.  This  is  in  agreement  with  the 
general  findings  of  other  researchers  (19,  20). 


3.0  TASK  2:  EXTENSIONS  TO  TRAPS  FOR  NOISE  ELIMINATION 


3.1  Introduction 


The  goal  of  Task  2  was  the  investigation  and  incor¬ 
poration  into  TRAPS  of  background  noise  statistics  relevant 
to  the  nuclear  monitoring  problem.  Since,  as  cited  in  Section 
2.1,  TRAPS  Stage  One  cancelled  all  signals  which  did  not  belong 
in  the  certain  class  of  transients  that  comprise  the  seismic 
events  of  interest,  initial  research  efforts  concentrated  on 
means  of  improving  TRAPS'  sophisticated  noise  elimination  stage. 
The  first  such  improvement  was  the  replacement  of  the  adaptive 
AR  model  of  the  original  TRAPS  with  an  adaptive  autoregressive 
moving  average  (AARMA)  model.  This  allowed  the  algorithm  to 
better  parameterize  both  the  signal  and  the  noise,  thus  in¬ 
creasing  TRAPS'  efficacy  in  a  low  signal  to  noise  (SNR)  environ¬ 
ment.  The  second  improvement  to  TRAPS'  Stage  One  was  the  ex¬ 
pansion  of  the  first  Stage  to  include  a  reference  waveform 
that  slaved  the  AARMA  model.  This  process  not  only  prevented 
jamming,  but  also  inhibited  the  passing  to  TRAPS'  Stage  Two 
noise  of  similar  characteristics  to  the  signals  of  concern. 

Both  these  improvements  to  TRAPS  were  detailed  in  the  semi¬ 
annual  report  and  shall  be  briefed  in  Section  3.2  for  the  sake 
of  continuity  in  this  final  report. 

The  research  efforts  in  the  first  six  months  produced 
an  accurate  AARMA  model,  which  was  constrained  to  remain  stable. 
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The  accuracy  and  consequential  noise  elimination  provided  by 
the  AARMA  model  for  the  majority  of  the  cases  is  impressive. 
Extensive  analysis  in  the  final  six  months  of  the  contract 
showed  a  weakness  in  the  model,  when  the  noise  had  a  sharp 
transition  in  its  characteristics.  While  this  weakness  was 
prevalent  for  only  a  short  time  span,  it  did  deteriorate  the 
model.  This  weakness  was  eliminated  by  putting  the  AARMA  into 
a  two  mode  process.  Mode  One  was  as  described  in  Section  3.2. 
Whenever  the  error  function  reflected  specific  conditions  the 
model  switched  to  Mode  Two  which  was  the  AARMA  mdoel  expanded 
into  a  cascade  structure  called  CAARMA.  Sections  3.3  and  3.4 
present  this  model  as  well  as  comparison  studies. 

Finally,  the  CAARMA  model  was  tailored  to  the 
nuclear  monitoring  problem.  The  noise  statistics  as  perceived 
by  an  array  were  evaluated  for  the  relationship  to  the  para¬ 
meters  of  a  CAARMA  model.  This  research  is  presented  in  Section 


3 . 2  An  Adaptive  Autoregressive  Moving  Average  (AARMA) 

for  Improved  Detection  in  Low  SNR  with  Anti-Jamming 
Strategies 


The  original  TRAPS  which  was  presented  in  Section 
2  was  based  on  an  adaptive  AR  process.  This  was  because  a 
priori  knowledge  about  the  generating  mechanism  of  the  time 
series  of  the  original  application  was  not  known.  However,  for 
the  situation  of  explosion/earthquake  monitoring,  the  generating 
mechanism,  as  well  as  the  background  noise  statistics  are  avail¬ 
able.  Cognizance  of  wave  propagation  models  along  with  initial 
evaluation  of  the  NORESS  data,  indicated  that  a  processing  gain 
would  be  achieved  by  upgrading  TRAPS  to  an  adaptive  ARMA  model. 
This  was  due  to  two  factors.  The  first  factor  is  that  the 
AARMA  generally  improves  spectral  resolution  at  low  SNR,  since 
the  ARMA  model  often  provides  better  spectral  estimates  than 
the  more  specialized  AR  or  MA  models  (21).  The  second  factor 
to  be  considered  in  applying  the  ARMA  model  to  seismic  signal 
processing  for  improved  detection  and  localization  is  that  as 
the  distance  from  source  to  receiver  increases,  there  is  the 
increase  in  multipath,  which  may  be  represented  by  an  ARMA 
system  (22).  The  ability  of  the  AARMA  to  improve  spectral  re¬ 
solution  and  to  accurately  model  the  NORESS  data  is  innately 
related  to  the  physics  of  the  problem,  which  we  shall  briefly 
outline  at  this  time. 


The  unprocessed  data  from  a  physical  point  of  view 
is  the  response  of  the  system  consisting  of  the  geophone  array 
and  the  media  through  which  the  impulsive  source  is  propagated. 
Thus  the  received  data  may  be  defined  as  the  sum  of  time- 
varying  responses  to  sharp  impulses,  each  sharp  impulse  asso¬ 
ciated  with  a  media  interface  that  reflects  the  energy  to  some 
degree.  The  time-varying  responses  are  the  wavelets.  The 
interpretation  of  the  received  waveform  consists  of  breaking 
down  this  elaborate  wavelet  into  its  component  wavelets.  The 
mathematical  solution  to  the  interpretation  problem  can  be 
understood  by  first  assuming  that  a  given  section  of  raw  data 
trace  is  additively  composed  of  wavelets.  This  wavelet  shape 
can  be  represented  by  a  stable,  casual,  minimum-delay  time 
function.  Thus,  each  wavelet  is  assumed  to  be  a  one-sided 
transient  that  damps  with  a  certain  degree  of  rapidity.  Also, 
knowledge  of  the  arrival  time  of  one  wavelet  does  not  allow 
prediction  of  the  arrival  time  of  another  wavelet.  It  is 
further  assumed  that  knowledge  of  the  strength  of  one  wavelet 
cannot  predict  the  strength  of  another  wavelet. 

These  seismic  trace  assumptions  can  now  be  trans¬ 
lated  into  mathematical  notation  for  discrete  time  t.  First, 
the  fundamental  constant  wavelet  can  be  represented  at  discrete, 
equally  spaced,  integer-valued  times  t.  Relative  to  the  nature 
of  this  wavelet,  its  shape  is  given  by  a  set  of  coefficients, 
bt,  that  are  weighted  by  their  respective  strengths,  it.  That 


is,  at  time  t,  the  wavelet  that  arrives  then  gives  the  con¬ 
tribution  i^b  ;  the  wavelet  that  arrived  at  t-1  gives  the 
t  o 

contribution  i  jb^;  etc*  The  seismic  trace,  xt,  at  time  t, 
is  then  the  summation  of  all  these  contributions.  This  wavelet 
complex  may  then  be  written  as  the  finite  sum 


xt  ■  Vt  *  Vt-l  4  b2it-2  * 


£  Vt- 

s*  0  s  1  • 
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for  the  arrival  times  between  t^  and  t£  called  the  time  inter¬ 
val  or  time  gate,  which  comprises  our  basic  section  of  seismic 
trace.  This  equation  includes  tails  of  wavelets  with  shape, 
bt,  these  wavelets  being  due  to  sharp,  yet  not  infinitely  sharp, 

impulses  i  -1,  i  -2,...  which  occur  before  time  t. .  Without 
Zl  Z1 

loss  of  generality,  the  sharp  impulses  it  may  be  centered  so 
that  their  mean  E{it\  is  equal  to  zero.  In  the  above  equation 
coefficients  b^  represent  the  "dynamics"  of  the  time  series, 
whereas  the  finite  duration  stochastic  function  it  represents  is 
the  "random"  nature  of  the  time  series. 

It  is  conceptually  easy  to  imagine  this  decomposition 
of  a  seismogram.  It  is  quite  another  matter  to  solve  the  prob¬ 
lem  of  separating  the  "dynamic"  and  the  "random"  components  of 
a  finite  section  of  data.  As  is  well  known,  state-variable 
methods  and  models  are  efficient  means  of  handling  time-varying, 
non-linear  and/or  nonstationary  effects,  multichannel  data, 
and  uncertainty  (23).  A  state-variable  model  is  that  set  of 


equations  which  describe  the  unique  relations  between  the 
input,  output  and  state.  The  state  variables  in  the  geophone 
and  layered  earth  system  has  a  well-defined  physical  meaning, 
as  just  described.  It  is  thus  possible  to  identify  the  system 
by  the  structure  of  the  transfer  function.  The  ARMA  model  is 
such  a  transfer  function.  An  ARMA  model  can  be  written  in  the 
following  notation: 

p  q 

x(n)  -£a.x(n-i)  =  XZ  b.W(n-i)  3.2 

i=l  1  i  =  0  1 

in  which  x(n)  is  the  modeled  time  series  and  W(n)  is  a  related 

"white"  noise  series.  Takinq  the  z  transform  (Z=exp  (i2irAtf )  , 
we  see  that  we  have 

BCZj 

X(Z)  =  JJV)  W(Z)  3.3 

♦ 

where  A(z)  and  B(z)  are  p,q  degree  polynomials  in  z,  and  where 
W(z)  is  the  transform  of  the  series  W(n) .  We  consider  the 
observed  series  x(n)  to  be  the  output  of  a  filter  with  transfer 
function  H(z)  =  B(z)A(z)  and  input  W(n). 

So  as  to  better  understand  the  system  identification 
by  the  structure  of  the  transfer  function,  the  transport  media 
may  be  viewed  as  made  up  of  layered  strata  and  the  seismic  wave 
propagation  in  this  layered  medium  as  reasonably  modeled  by  plane 
waves.  The  matchinq  of  wave  equation  jump  conditions  at  each 
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interface  leads  to  the  definition  of  a  reflection  coefficient, 
as  discussed  earlier  in  this  section.  Thus,  the  response  of 
a  system  of  varied  density  to  a  unit  spike  is  given  by  the 
equation  3.2,  where  the  a's  are  related  to  the  reflection  co¬ 
efficients  and  the  b's  are  functions  of  the  media  associated 
with  factors  that  could  be  considered  as  interfering  with  an 
accurate  event  detection.  We  will  consequently  categorize  the 
variations  in  the  signal,  which  contribute  to  these  terms,  as 
noise.  The  ARMA  equations  need  now  to  be  specified  over  a 
finite  time  and  devise  a  method  for  solving  the  unknowns  in 
an  environment  which  corrects  for  gradual  variation  of  the 
transient  signal.  Referring  to  our  earlier  discussion  of  the 
advantages  provided  by  an  adaptive  process  in  Section  2,  we 
propose  to  use  an  adaptive  ARMA  model.  If  we  were  concerned 
only  with  an  adaptive  AR  process,  the  equations  to  be  solved 
would  be  quadratic  in  nature  and  the  unknowns  could  easily  be 
resolved,  as  discussed  in  Section  2.1.  The  adaptive  ARMA 
model  is  nonlinear  and  thus  requires  a  method  which  remains 
stable  within  the  process  of  solving  the  nonlinear  equation 
without  requiring  excessive  computer  time.  In  order  to 
accomplish  this,  the  following  hybrid  adaptive  ARMA  model  was 
developed.  First,  the  model  in  equation  3.2  is  rewritten  in 
two  stages  as  follows: 


I:  e(n)  H  x(n)-$I  a.x(n-i) 

i  =  l  1 
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<« 


0 


and 

P 

II:  7"  b.W(n-i)  =  e(n)  3.5 

i-0  1 

It  is  possible  to  initially  calculate  an  estimate, 
(n) ,  for  the  series  e(n)  (2).  This  leaves  the  problem  of 
estimating  b^,  Wn  given  equation  3.5.  We  now  estimate  b^  and 

A  A 

Wfn)  as  those  values  b^,  W(n)  which  minimize  the  total  squared 
error : 

Gn(b,W)  =IlFn(b,W)-a(n)^2  3.6 

Subject  to  the  constraints 

_  C  p=0 

m  r 

R  (n)  =  YL  W(n- i)W(n-p-i)  =■ 

P  i=n 

0  l^p^scm  3 . 7 

where  C  is  an  arbitrary  constraint,  s  is  the  number  of  unknown 
W  and 

A/ 

q 

F  (b,W)  =  £  b.W(n-i)  3.8 

n  i=0  1 

Note  that: 

1)  At  this  stage  &  (n)  is  considered  given. 

2)  Any  algorithm  obtaining  £(n)  must  eventually 
be  justified 
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3)  More  precision  on  limits  to  sums  will  be 
made  subsequently. 

We  now  recast  these  equations  into  a  form  suitable  for 
numerical  treatment. 

Suppose  that  at  time  tn*nAt  we  have  observed  (estimate) 
e(n).  (We  drop  the  carets  for  notational  convenience.)  We 
compute  F  from  a  backward  sum 


F_.k  =  Z  b (i+l)W(n-k- i) 
n  i-0 


where  k=0,...,N  (n>N) 

and  b(i  +  l)  =  b^  where 


W(n-k-i)  =  Wnk_.  (n>k+i) 
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are  FORTRAN  arrays  representing  the  unknowns.  We  represent 
these  unknowns  jointly  as  f  (at  epoch  n) : 


lW (n-m+1) 


m  <  s 
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lb  (m-s)  m  >  s 


where  l<m<s+q+l  and  q+1  is  the  number  of  unknown  b. 

—  —  14  A* 

We  write  our  problem  (at  epoch  n)  as 


ML  (Fn-k<?>-vi/  *  LVVCV  ■ 0  3-12 


■)- 
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where  A  are  Lagrange  multipliers. 

P 

The  normal  equations  are: 


2B  K-x-v*5 


3F  .  s-1  5Rp 

+  r  L  - 1-  0  3.13 


9$m  p*0  p  3fm 


Rp  (fg)  =  C6 
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These  equations  must  be  solved  for  £  (using  the  constraints 
to  be  X  ),  but  are  nonlinear.  We  use  the  Gauss-Newton  itera¬ 
tion  method  which  ignores  certain  second  derivatives  to  obtain 


zb  ]  Cf1*1-,1)* 

t  k=0  9£m  <^jlt  p=0  p  ^fm^^t  J  ^ 


s-1 

£  r  li+ili 


P=0  ^m 


-V  '  *  [  2  £.  'WW' 


n-k  + 


s-1 

£  X, 
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where  i  refers  to  the  iteration  number. 

And  the  constraints  are  linearized  to 


r  (pi-1  -  r  )  =  -  1 
a?  t  t  ?  t J  Rp 


3.16 
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These  equations  may  be  written  in  matrix  form  as: 


Ea^=I 

A.  J-  - - '  ~ 

^  A- 


3.17 


HAf  =  J 


3.18 


where 


er  .  2  ?  3Iz±  i!nri  t  siS  J*!be 

mt  k=o  3?m  3$t  p=o  Ap 


mp 


dRp 


m 


"I 


N 

2  £  (Fn-k 

k=0  n  k 


e  ,  ) 
n-k' 


j-JL-jF.  +  V;1  1  & 

m  p=0  p  3?t 


h  -  i^E 
HPt  -  3it 


J  =  -Rp 
o 


A  single  matrix  equation  may  be  written  as 


O' 

g 

s 

'  £  ' 

H 

0J 

*1 

.  1  . 

3.19 
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This  matrix  equation  is  solved.  And  the  results  are  used  to 
update  the  iterates 

fn+1  -  fn  +a£  .  etc.  3.20 

Despite  the  excellent  detection/classification 
capability  of  the  original  TRAPS,  due  to  its  adaptive  nature, 
it  did  exhibit  signal  cancellation  phenomena  like  any  adaptive 
filter  (24-26,17),  when  adapting  rapidly.  This  phenomena  is  a 
consequence  of  the  adaptive  interaction  between  the  desired 
signal  and  interference.  This  interference  could  be  anything 
from  bandpass  noise  to  a  sum  of  sinusoids  suitably  spaced  in 
frequency.  The  existence  of  these  problems  in  the  AARMA  model 
as  well  as  the  AR  model  motivated  the  development  of  a  more 
complex  first  stage  of  TRAPS.  The  method  developed  has  simili¬ 
tude  with  work  of  Ken  Duvall  (27). 

A  schematic  of  the  adaptive  noise  canceller  at  its 
development  stage  per  the  semi-annual  report  is  shown  in  Figure 
8.  The  "primary  input"  contains  a  useful  signal  s,  plus 
interference  nQ.  The  "reference  input”  contains  interference 
n^,  related  to  that  of  the  primary  input,  and  is  separately 
obtained  in  practical  systems.  The  relationship  between  the 
two  interferences  is  generally  unknown  a  priori .  The  adaptive 
filter  has  the  job  of  shaping  the  reference  interference  to 
replicate  (in  the  least  squared  error  sense)  the  primary 


interference  so  that  subtraction  will  remove  the  interf erencre 
from  the  primary  input  and  thereby  deliver  a  much  mere  useful 
output.  It  has  been  shown  in  (28)  that  an  adaptive  filter 
minimizing  output  power  in  the  system  of  Figure  8  causes  the 
system  output  to  be  a  best  least  squares  estimate  of  the  use¬ 
ful  signal  s.  Unfortunately,  under  certain  circumstances , 
primary  input  components  may  be  partially  cancelled  and  dis¬ 
torted  despite  the  fact  that  they  are  not  correlated  with  the 
reference  input.  This  distortion  in  many  cases  would  be 
quite  acceptable.  What  is  not  acceptable  are  the  "non-Weiner" 
effects  due  to  weight  fluctuation  which  give  rise  to  signal 
cancellation. 

In  order  to  facilitate  discussion  of  the  solution  to 
the  problem  of  these  non-Weiner  effects,  consider  the  example 
where  a  sinusoid  jammer  from  a  direction  other  than  that  of 
the  primary  signal,  or  "off  the  look  direction",  causes  sinu¬ 
soidal  fluctuations  in  all  the  weights.  The  look  direction 
constraint  is  unit  gain  and  zero  phase  over  all  frequencies. 
The  useful  signal  arriving  in  the  look  direction  encounters 
a  unit  gain  with  the  main  beam  so  constrained.  The  sinusoid 
off  the  lock  direction  for  this  example  produces  a  notch 
along  the  primary  signal  flow  path  to  the  output  by  means  off 
the  fixed  weight  filter.  Notching  phenomena  in  this  system 
are  very  much  like  those  of  the  adaptive  canceller  of  Fig.  8. 
Signals  from  the  look  direction  do  not  appear  at  the  adaptive 


filter  inputs;  only  interference  is  present  there.  Both 
signal  and  jammer  are  present  in  the  primary  signal  flow  path, 
and  both  signal  and  jammer  experience  notching  at  the  jammer 
frequency.  With  high  speed  adaptation,  there  is  the  danger 
of  eliminating  the  signal  and  jammer.  The  signal  cancellation 
phenomena  that  has  just  been  discussed  is  due  to  interaction 
between  the  signal  and  the  jammer  in  the  adaptive  beamformer. 
Since  the  interaction  is  the  root  of  the  problem,  it  is  useful 
to  consider  adaptive  structures  that  separate  the  signal  and 
the  jammer,  during  the  adaptive  process.  A  remedy  to  the 
signal  cancellation  problem  will  now  be  discussed.  The  method 
devised  is  based  on  original  work  by  Ken  Duvall  (27) .  The 
Duvall  system  incorporates  two  signal  processing  systems,  one 
to  perform  the  adaptation,  the  other  to  generate  the  system 
output  signal. 

The  system  incorporated  into  TRAPS'  Stage  One  is 
shown  in  Figure  9.  It  makes  use  of  two  filters  per  phone. 

The  filter  on  the  right  is  connected  directly  to  the  elements 
and  is  used  to  derive  the  array  output  signal.  It  is,  however, 
a  slaved  filter  rather  than  the  adaptive  filter  that  would 
usually  be  expected  in  this  position.  The  filter  on  the  left 
is  the  TRAPS  adaptive  filter  that  is  connected  to  the  array 
elements  through  a  subtractive  preprocessor.  The  preprocessor 
excludes  the  look-direction  signal,  but  admits  ja"’  °r  signals 
in  a  modified  form.  The  adaptive  process  generates  a  set  of 


weights  that  provides  some  specified  look-direction  gain  while 
minimizing  (in  the  least-squares  sense)  the  jammer  contribu¬ 
tion.  These  weights  are  then  input  into  the  slaved  beamformer 
to  provide  the  desired  signal  reception  and  jammer  rejection. 

This  study  of  a  slaved  first  stage  provided  insight 
into  a  Stage  Two  system  which  not  only  prevented  the  signal 
cancellation,  but  gave  the  algorithm  perspective  as  to  the  path 
the  coefficients  were  to  take,  once  the  error  exceeded  a 
specific  threshold.  This  system  for  Stage  Two  was  labeled 
"cascade  structure"  and  is  detailed  in  Section  3.3. 


3 . 3  Stage  Two  Cascade  Formalism 

What  the  algorithm  detailed  in  Section  3.2  permits 
is  the  obtainment  of  an  accurate  ARMA  model,  which  is  con¬ 
strained  to  remain  stable.  The  accuracy  of  the  ARMA  model 
for  the  majority  of  the  cases  is  impressive.  Extensive  analy¬ 
sis  since  the  semi-annual  report  showed  a  weakness  in  Stage  Two 
of  TRAPS,  when  the  noise  or  signal  had  a  sharp  transition  in 
its  characteristics.  While  this  weakness  was  prevalent  for  only 
a  short  span  of  time  and  did  not  deteriorate  the  model  at  high 
SNR,  it  did  reduce  localization  accuracy  at  low  SNR.  This 
problem  was  therefore  resolved  by  what  we  shall  term  a  "cascade 
structure."  The  primary  goal  of  this  section  is  to  detail  the 
development  of  a  cascade  structure  to  the  AARMA  model  of  TRAPS’ 
Stage  Two.  So  that  the  processing  gains  allowed  by  this  cas¬ 
cade  structure  can  be  best  understood,  this  section  will  also 
brief  the  characteristics  of  the  mean-squared-error  surface 
and  the  performance  of  a  recursive  LMS  algorithm  for  use  in  the 
evaluation . 

A  performance  surface  in  terms  of  the  mean-squared- 
output  error  is 

E(e2)  =  E[(y-y)2]  =  E(y2)  +  E  (y2)  -  2E(yy)  3.21 


So  that  the  association  between  the  cascade  structure  and  the 


improvements  to  the  AARMA  may  be  easily  understood,  the 
developmental  steps  will  be  restricted  to  the  case  where  the 
model  and  the  adaptive  algorithm  are  both  second  order  with 


H(z) 


Btzl  .  bQ  *  V'*  +  b2* 

1  -  a1z"1  -a2z-2 
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The  assumption  that  y  and  y  are  stationary  random  signals  with 
zero  mean  and  finite  variance,  and  that  the  input  u  is  white 
noise  with  unit  power  density,  allows  a  normalized  error  surface 
to  be  defined  in  terms  of  the  transfer  functions  [29] 


1  + 


2trj«- 


“7  §  H(z)  [H*  (z)  -  2H*  (z)  ]~— 
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where 

*  ’  WJ  j  [H(z)]2 
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A 

H(z)  is  the  transfer  function  of  the  adaptive  algorithm,  and 
H*(z)  denotes  the  complex  conjugate  of  H(z). 

Attention  shall  now  be  confined  to  the  two-pole  model 
in  order  to  illustrate  error  surface  effects  with 


H(z) 


l-afz  ^ 
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This  allows  the  mapping  of  the  interior  of  the  unit  circle  in 
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the  z-plane  onto  a  triangular  region  in  the  plane  associated 
with  the  two  adaptive  coefficients.  The  interior  of  the  tri¬ 
angle  shown  in  Figure  10  thus  represents  the  stable  region  for 
the  adaptive  coefficients  a^, 

Plots  illustrating  the  contours  of  constant  mean- 
squared-error  for  three  different  models  are  illustrated  in 
Figure  11.  These  plots  emphasize  the  extreme  nonuniformity  of 
the  ARMA  performance  surface.  The  contours  tend  to  be  elongated 
in  shape,  so  that  the  gradient  vector  does  not  always  direct 
itself  toward  the  minimum  where  £  ■  0 .  A  more  complete  analysis 
of  the  error  surface  of  AARMA  algorithms  is  provided  in  reference 
[29]. 

The  parameter  update  formula  for  an  adaptive  co¬ 
efficient  vector  9  is  expressed  as 

Gk+1  =  0k  +  A (”vk^  3,26 

where  jm  is  a  parameter  which  controls  stability  and  rate  of 
convergence,  and  is  the  value  (or  an  estimate)  of  the  error 
surface  gradient  where  0  =  0^.  Detail  of  this  can  be  found  in 

Section  2.1. 

For  the  comparisons  to  follow  it  is  important  to 
realize  that  for  a  fixed  point  on  the  mean-squared-error  sur¬ 
face,  the  estimate  obtained  from  the  recursive  LMS  algorithm 
is  unbiased.  To  see  this,  we  define  a  coefficient  vector 


A.  A  A  A  A  A  j 

0  =  [ax  a2  b0  bx  b2]1  . 


3.27 


The  true  gradient  vector  is  therefore  expressed  as 


k  3  0 


3.  28 


where 


?k  ■ 
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The  gradient  estimate  for  the  recursive  LMS  algorithm  is  derived 
in  terms  of  the  instantaneous  squared  output  error  as 


A 

V, 


3£k 

,  A 
^0 


3.30 


where 


3.31 


We  now  define  the  expected  value  of  the  gradient  estimate  over 
k  as  the  average  of  an  arbitrarily  large  number  of  samples: 
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3 . 32d 
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3 . 32e 

• 

The  fact  that  ECv^)  =  v  k  implies  that  the  gradient  estimate  is 
asymptotoically  unbiased. 

Simulation  Results: 

Again  we  will  restrict  our  attention  to  the  two-pole 

A* 

model  and  fix  the  numerator  coefficients:  =  bQ  *  1 ,  » 

/\  A 

b2  =  bi  =  b2  =  0.  In  Figure  32,  parameter  track  plots  illus¬ 
trate  the  performance  of  the  recursive  LMS  algorithm  for  each 
of  three  error  surfaces. 

It  is  apparent  that  the  gradient  estimates  are  very 
noisy  in  these  examples.  In  a  nonideal  environment,  where 
additional  noise  is  introduced  by  signal  measurement  and  other 
uncontrollable  sources,  generally  this  amount  of  estimate  noise 
would  be  unacceptable.  Smoother  parameter  tracks  may  be  ob¬ 
tained  at  the  expense  of  increased  convergence  time  by  reducing 
the  size  of  n.  This  is  illustrated  in  Figure  13.  The  recur¬ 
sive  LMS  algorithm  did  converge  in  every  example  when  p  was 
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properly  selected.  When  some  a  priori  knowledge  of  the  model 
characteristics  is  not  available,  choice  of  an  appropriate 
value  remains  an  empirical  process. 

From  the  preceeding  examples,  we  see  that  a  constant 
value  of  n  does  not  work  well  for  the  extremely  nonuniform 
error  surfaces.  Ideally,  we  would  like  to  take  small  steps 
when  the  gradient  is  very  steep  to  avoid  bouncing  from  side  to 
side  across  the  "valley."  If  n  is  so  selected,  when  the  co¬ 
efficients  reach  the  floor  of  the  valley,  convergence  to  the 
minimum  will  be  very  sluggish  due  to  the  small  gradient  [29]. 

The  definition  of  a  cascade  structure  can  now  be 
attempted.  The  basic  procedure  for  Stage  Two  will  still  follow 
the  technique  described  in  Section  3.2.  When  e  reaches  a  pre¬ 
determined  upper  limit,  the  algorithm  switches  to  a  cascade 
structure.  When  the  cascade  structure  is  in  position  the 
transfer  functions  are  of  the  form 

N 

H(z)  =  TT  Hj(z)  3.33 

where 

3.34 

3.35 
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The  adaptive  cascade  structure  is  illustrated  in  Figure  14. 

Applying  the  steepest  descent  criteria  defined 
earlier,  the  coefficient  updates  are  obtained  of  the  form 


Vk+’> 


b^j(k)  +  y(- 


a^jtk+l)  -  a4<(k)+r\|- 


1J 


*LY 

II 

2 

j 

CNJ 

m 

11 

SbJ 

1  =  0, 

1,  2 

j 

»  N 

1=1, 

2 

j 

■1,2 

3au  / 

.  N-l 


3 . 36a 


3.36b 


As  in  the  derivation  of  the  recursive  LMS  algorithm, 
we  now  assume  that  the  instantaneous  output  error  is  a  local 
estimate  for  the  mean-squared  error  to  obtain 

<$-  =  2e(k,  i|l£l 

3bij  9bij  3.37a 


2e(k) 


3a 


1j 


3.37b* 


Now  we  make  use  of  the  fact  that  the  model  output  y  is  inde¬ 
pendent  of  the  adaptive  coefficients  and  define  the  partial 
derivatives 


efj(k) 


3.38a 
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ij 


Combining  Equations  3.36  -  3.38,  the  coefficient  updates  are 
obtained  as  follows 


bi  j (k+1 )  3  ^j(k)  +  ue(k)g^(k) 


1  *  2  j  =  1,  2  ...  N-l 

j  =  0.  1,  2  j  =  N  3.39a 


a1j(k+1)  3  ajjOO  +ne(k)aj  (k)  i  =  i,  2 


J  3  1,  2  ...  N. 


3.39b 


In  order  to  derive  expressions  for  and  flN 

ij  ^  13 

defined  in  Equation  3.38,  we  will  make  use  of  a  series  of 
intermediate  values 


B"  (k)  4  ^ 
3bij 

a  3Vk) 


a1 J (k) 
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n  =  1  ,  2  ...  N 


n  =  1 ,  2  ...  N 


3.40a 


3.40b 


1j 


Since  the  output  of  stages  1,  2  ...  j  -  l  is  unaffected  by  the 
coefficients  of  stage  j,  we  see  that 


B"j(k)  -0  n  <  j 


0  n  <  j 


3.41a 


3.41b 
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When  n  =  j  ,  the  values  are  obtained  exactly  as  for  the  second 
order  section  in  Section  2. 1,3. 3  with  the  input  now  equal  to 

(or  u(k)  if  j=l)  and  output  .  Suppressing  the  time  de¬ 
pendence  of  the  adaptive  coefficients  for  notational  con¬ 
venience  we  obtain 


BjjOO  •  ♦^J<*-Z) 


°1j<k>  *  »'j{k-,)  *  *  "2jatj(k-2) 


3.42b 


For  n  >  j ,  the  partial  derivatives  are  formulated  as  filtered 
versions  of  p  ^  and  j .  : 


6U<k>  ■  V!j'<k>*iln^'lk-'>‘i2n8j],U-2) 
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+  alna1j(k_1 J  +  a2naij(k_2) 


3.43b 


Note  that 


n  <  N 


3.43c 
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If  we  now  assume  that  small  changes  are  made  at 

A  A 

each  iterative  cycle  so  that  we  obtain  a 

slightly  simpler  version 

i  =  1,  2  j  =  1,  2  ...  N-l 

i  =  0,  1,  2  j  =  N 

3.44a 

®ij  (k+l)  =  j  (k)  +  -»v.e  (k)«^  (k- i)  i  =  1,  2  j  ■  1,  2  ...  N-l 

3.44b 

N  N 

where  and  e<  j  are  formulated  as  illustrated  in  Figure  15.  An 
example  illustrating  the  convergence  of  a  four  pole  system  is 
provided  in  Figure  16. 


N 


bi;j(k+l)  =  b.^k)  +  /x  e  Ck^/3  j  (k-  i) 


3 . 4  Evaluation  of  Cascaded  AARMA  on  Simple  Simulations 


A  small  class  of  AARMA  algorithms  with  comparable 
computational  complexity  and  proven  convergence  characteristics 
will  be  used  to  demonstrate  the  capability  of  the  cascade 
structure.  By  choosing  the  order  of  the  adaptive  algorithm  to 
be  at  least  as  high  as  that  of  the  model,  a  unimodal  perfor¬ 
mance  surface  is  assured  [29], 

In  addition  to  the  adaptive  algorithms  presented 
earlier,  we  will  employ  the  SHARF  algorithm.  The  analysis  by 
Ljung  [30-31]  provides  a  convergence  proof  for  the  RLMS  algo¬ 
rithm.  Since  SHARF  is  not  a  gradient  search  algorithm,  the 
convergence  proof  relies  upon  the  concept  of  hyperstability 
[32-36] . 

A  direct  form  realization  of  the  transfer  function 
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is  utilized  by  both  the  RLMS  and  SHARF  algorithms. 

Since  the  second  order  RLMS  algorithm  discussed  in 
the  previous  Section  was  constrained  to  remain  stable  during 
the  adaptive  process,  use  will  also  be  made  of  the  general  Nth 
order  unconstrained  realization.  This  then  provides  a  wide 
range  of  algorithms  for  the  following  examples.  The  adaptive 
updates  are  given  by 
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b^(k+l)  =  b^(k)  +  m  e  (k)y?  (k- i)  i  =  0,  1  ...  N  3.46a 


a^(k+l)  =  a^(k)  +  xe  (k)*.(k-i)  i  =  1,  2  ...  N  3.46b 


where  oC.yS  and  e  are  formulated  as  in  Figure  15. 

Use  of  the  SHARF  algorithm  requires  the  selection  of 
a  set  of  error  smoothing  coefficients,  c^,  ...  Cp  [32].  In 

order  to  satisfy  the  hyperstability  criterion,  these  parameters 
must  be  chosen  so  that 


3.47 


where 


C(z)  =  1  +  c1z’1+CpZ  p 


3.48 


and  A(z)  is  the  denominator  of  the  fixed  transfer  function  H(z)'. 
As  a  result,  some  a  priori  knowledge  of  the  model  is  required. 

An  alternate  technique  has  been  proposed  which  allows  recur¬ 
sive  estimation  of  the  C(z)  coefficients  [37].  Here  we  will 
use  the  fixed  C(z)  polynomial.  The  adaptive  updates  for  the 
SHARF  algorithm  may  be  expressed  as 


b^(k+l)  =  b^(k)  +u.u(k-i)  v(k)  i  ■  0,  1  ...  N  3.49a 
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t 


i  ■  1 ,  2  .  . .  N  3 . 


ai(k+l)  =  ai(k)  +  *.y(k-i)  v(k) 


where  u,  y  and  v  are  formulated  as  in  Figure  15  [33]. 

A  second  order  model  with  transfer  function 


H  (z)  = 


.  2871 


4466z 


3.50 


1.  -  1.2971a  +  . 6949z 


will  be  the  first  example  to  be  used  for  comparison.  The  plots 
in  Figures  17  and  18  illustrate  the  reduction  in  output  error 
versus  the  number  of  data  samples  used  for  four  different 


cases . 


In  Figure  17. a,  the  RLMS  algorithm  defined  by  3.46 


was  used.  With  convergence  factors »  .07,  the  resultant 
steady-state  transfer  function  was 


H(z) 


RLMS 


.2871  -  .  4466z  -  .OOOlz 

1.  -  1 . 2970z" 1  +  .6949z"2 


3.51 


Figure  17. b  illustrates  the  performance  of  Cascade 
AARMA  (CAARMA)  algorithm  defined  by  (27).  Using  convergence 
factors  *  .5,  the  steady-state  transfer  function  was 


H(z) 


CAARMA 


.  2871  -  .  4466z  x 
1.  -  1.2971z_1  +  . 6949z 
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In  Figure  18,  the  SHARF  algorithm  defined  in  (3.49) 
was  used.  For  the  transfer  function  defined  by  (3.50),  we 
can  show  that 


Re  l  sferl  > 
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Thus,  C(z)  *  1  so  that  no  error  smoothing  is  performed. 
Figure  18. a  illustrates  the  performance  of  SHARF  for  this 
case  with  convergence  factors  =  JjL  =  .07.  The  steady-state 
transfer  function  was 


H(z)SHARF 


*  -2871 
1.  - 


-  .4466z 
1.2970Z-1  + 


■  OOOlz 
.  6947z~ 


3.54 


It  has  been  shown  that  choice  of  the  C(z)  coefficients  also 
affects  the  rate  of  parameter  convergence  [32].  This  is 
illustrated  by  Figure  18. b.  Here  we  let 

C(z)  =  1.  -  1 . 25z_1  3.55 


and  again  use  convergence  f actors  7^  =,^C  =  .07.  The  resultant  » 
steady-state  transfer  function  was 


H(z^SHARF 


.2871  -  . 4466z-1 
1.  -  1 . 2971z  1  +  . 6949z 


3.56 


Note  the  improvement  in  the  rate  of  parameter  convergence  which 
results  from  the  error  smoothing.  For  further  analysis  of  the 
relationship  between  the  C(z)  polynomial  and  the  rate  of  con¬ 
vergence,  see  References  [32-36]. 
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From  the  results  of  these  simulations,  we  see  that 
the  CAARMA  algorithm  evidenced  the  best  performance  in  terms 
of  the  number  of  data  samples  required  for  coefficient  con¬ 
vergence.  An  attempt  was  made  to  maximize  the  rate  of  con¬ 
vergence  in  all  cases  by  adjusting  -^and^ix,  but  as  was  remarked 
previously,  SHARF  convergence  characteristics  may  be  improved 
by  a  different  selection  of  C(z). 

It  has  been  shown  that  when  the  order  of  the  adaptive 
algorithm  exceeds  that  of  the  model  and  the  performance  surface 
is  unimodal ,  the  minima  is  distributed  [29],  A  fourth  order 
adaptive  algorithm  will  now  be  used  as  an  example  in  an 
attempt  to  identify  a  third  order  model  with  transfer  function 


H  (z) 


.0565(1.  +  z‘1)(l.  -  1.0166Z"1  +  z~2) 

(1.  -  .6830z  ■1")(1.  -  1.4461Z-1  +  .7957z"2) 


This  is  the  digital  equivalent  of  a  low-pass  elliptic  filter 
which  was  obtained  via  the  bilinear  transformation  [38].  The 
plots  in  Figures  19  and  20  illustrate  the  reduction  in  output 
error  versus  the  number  of  data  samples  for  four  different 
adaptive  algorithms. 

The  cascade  structure  algorithms  were  used  in 
Figure  19.  Figure  19. a  illustrates  the  performance  of  the 
CAARMA  algorithm  defined  by  (3.42).  With  convergence  factors 
77  * .025,  the  steady-state  transfer  function  was 


H(z) 


CAARMA 


f  1-  ~  1.0164Z'1  +  .  9999z  ~  2  \ 

1.  -  1.4462Z'1  ♦  .  7  9  5  8  z  "  2 
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3.58 


The  second  stage  can  be  factored  as 


,0562(1.  ♦  . 998z  ) (1  +  .370z 

(1  -  .683z_1)  (1  +  . 371z -1) 


3.59 


thus  illustrating  the  existence  of  a  common  factor  in  H(z) 
which  results  from  the  distributed  minima. 

The  two  direct  form  adaptive  algorithms  defined  at 
the  beginning  of  this  Section  were  employed  during  the  simula¬ 
tion  illustrated  in  Figure  20.  The  direct  form  equivalent  of 
the  model  defined  by  (3.57)  is 


H(z) 


.0563(1.  -  . 0166z 


0166z"2  +  z'3: 


2.1291Z  1  +  1.7834z  ‘  -  .5435z 


3.61T 


Figure  20. a  was  obtained  via  the  RLMS  algorithm  defined  by 
(3.46)  with  convergence  factors  .Ol,^*  .005.  The 
resultant  steady-state  transfer  function  can  be  expressed  as 


H(z) 


RLMS 


.0562(1.  ♦  .867z'1)(l.  -  .0165Z'1 


. 0165z  * 2  +  1 . 002z"3 


(1.  +  .867z"j  (1 .  -  2 . 130z " 1  +  1.7836Z  -  .5432z*°) 


3.61 
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thus  demonstrating  the  existence  of  a  common  factor  due  to  the 
distrubuted  minima.  A  fourth  order  SHARF  algorithm  was  em¬ 
ployed  during  the  simulation  illustrated  in  Figure  20. b.  To 
ensure  satisfaction  of  the  strictly  positive  real  condition, 
we  let 


C(z)  -  1.  -  1.4461Z-1  +  . 7957z 


3.62 


to  obtain 


R6^]  '  ReL~-,6830z-^>  “  ' 

With  convergence  factors  ju.  »  .07,  the  resultant  steady-  * 

state  transfer  function  was 


H^z)SHARF 


.062(1.  +  .  762z~ 1-)  (1 .  -  .  0164Z1  -  ,0165z~2  ♦  .9998z~5-) 
(1.  +  .762z'1)(l.  -  2 . 1 2 8z" 1  +  1.7836z'2  -  .5434z"^) 


3.64 


This  example  serves  a  dual  purpose.  It  illustrates 
the  performance  of  the  multistage  cascade  algorithms  and  demon¬ 
strates  that  when  the  performance  surface  minima  is  distributed, 

A 

the  transfer  function  H(z)  contains  a  common  factor. 

It  is  evident  from  the  illustrations  in  Figures  19 
and  20  that  the  cascade  structure  algorithms  yielded  the 
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superior  rate  of  convergence  under  the  conditions  previously 
described.  Maintaining  stability  during  the  adaptive  process 
was  found  to  be  the  primary  limiting  factor  for  both  RLMS  and 
SHARF.  It  is  possible  that  a  different  selection  of  C(z) 
could  improve  the  rate  of  convergence  of  SHARF,  however  the 
choice  given  by  (3.62)  already  assumes  a  gr^at  deal  of  a 
priori  knowledge  of  model  parameters. 

An  idealized  operating  environment  has  been  used  in 
the  previous  simulations,  so  as  to  yield  insight  into  the 
characteristics  of  this  class  of  adaptive  algorithms.  The 
performance  of  the  adaptive  algorithms  presented  earlier  will 
row  be  examined  under  nonideal  operating  conditions. 

The  first  condition  that  will  be  examined  is  that 
of  the  algorithm  response  to  insufficient  model  order.  Guaran¬ 
teed  convergence  to  the  global  minima  requires  that  the  per¬ 
formance  surface  be  unimodal  when  gradient  search  techniques 
are  employed.  It  has  been  shown  that  when  the  order  of  the 
adaptive  algorithm  is  less  than  that  of  the  model,  the  perfor¬ 
mance  surface  may  contain  multiple  minima  even  though  the 
system  identification  configuration  is  utilized  and  the  input 
is  white  noise  [29].  In  this  case,  steady-state  adaptive  co¬ 
efficients  are  dependent  upon  their  initial  values  [30]. 

Since  the  performance  surface  is  quadratic  with 
respect  to  the  numerator  coefficients,  pole  adaptation  is  our 
primary  concern.  Here  we  will  assume  a  sixth  order  all -pole 
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model  with  transfer  function 


H(z)  °  _ _  .1 _ _ _ 

(1.  -  1.2z_1  7  .6z'Z)(l.  +  1-2Z'1  +  .6z"Z)Cl.  -  .64z'2) 

3.65 

The  corresponding  pole  locations  are  symmetrically  distributed 
in  the  z-plane  as  illustrated  in  Figure  21. 

Simulations  were  performed  using  a  fourth  order 
CAARMA  algorithm  with  three  different  sets  of  initial  conditions 
on  the  stage  two  coefficients.  Zero  initial  conditions  were 
used  for  stage  one  coefficients.  The  resultant  steady-state 
pole  locations  are  illustrated  in  Figure  22.  Plots  demonstrat¬ 
ing  the  reduction  in  output  error  for  each  case  are  provided 
in  Figure  23.  Note  that  the  error  reduction  in  Figure  23. c 
was  minimal,  thus  clearly  illustrating  the  primary  hazard  of 
multimodal  performance  surfaces. 

9 

The  next  example  that  shall  be  considered  is  that 
of  delay  adjustment.  Earlier  it  was  pointed  out  that  even 
when  the  order  of  the  adaptive  algorithm  is  sufficient,  the 
cascade  structure  defined  by  Figure  14  may  not  allow  the 

A 

unique  solution  H(z)  =  H(z)  because  of  a  delay  mismatch.  By 
using  the  alternate  structure  illustrated  in  Figure  24, 
however,  this  problem  can  be  alleviated. 

Consider  a  fourth  order  model  with  transfer  function 
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H(z) 
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Utilizing  a  fourth  order  cascade  structure,  we  see  that 


A  -1  A  -2 


B  (z)  =  Cl  +  bnz  +  b^z  A)(b02  +  b12z  +  b 2 2 z  *) 


3.67 


A  A 


A  A 


A  A 


=  b02  +  ^bllb02  +  b12^Z  +  (b21b02  +  bllb12 

A  A  A  AA  aA 

+  b22)z  +  (b21b12  +  b11b22)z  +  (b2ib22^z 

3.68 

A 

where  at  lease  one  of  the  stage  two  b  coefficients  must  be  non- 

A 

zero.  The  B(z)  polynomial  defined  by  (3.68)  cannot  match  the 
numerator  of  the  transfer  function  in  (3.66).  If  we  now  in¬ 
corporate  J.  units  of  pre-delay  as  illustrated  in  Figure  24, 
we  obtain  the  adaptive  numerator  polynomial 


_  0  A  A  -J?  A  A  A  A  A  A 

z  B (z)  *  bQ2z  +  (b11b02  +  b12)z  +  (b2ib02  +  bllbll  + 


A.  A  AA  AA 

u  ^  u  u  \  ^  .  /■  u 


b22-)z  2  +  (-b21b12  +  bllb22-)z 


+  Cb21b22)z 


-X- A 


3.69 


From  (3.69)  we  see  that  we  can  now  match  the  model  B(z)  in 
three  cases: 


z" 2B (z)  -  B (z)  ;  b. 


3.70a 


-  3  A  A 

z  b(z)  «*  B(z)  ;  b. 


3.70b 
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.5 


-4A 

z  ^B(z) 


=  B  (z) 


02 


3.70c 


A 

where  all  unspecified  coefficients  are  zero. 

Simulations  were  performed  for  two  cases  using  the 
fourth  order  CAARMA  algorithm  with  >*•  =  =  .5.  Figure  25 

illustrates  the  results  obtained  with  four  units  of  pre-delay. 
The  steady  state  transfer  function  was 


H(z) 


(1.  -  .  003z~ 1  +  .OOOz'2)  (.500  +  .  0041  z '  ^  -  .004z~2) 
(1.  -  1.199Z1  +  .600z  Z)(l.  -  1 . 008z  1  +  .401z  Z) 


3.71 

The  resultant  numerator  polynomial  can  be  expressed  as 

z‘4B(z)  =  . 500z"4  +  .002z"5  -  .004z'6  .  3.72 


Two  units  of  pre-delay  were  used  during  the  simulation  illus-  • 

trated  in  Figure  25. b.  The  steady-state  transfer  function  was 


H(z) 


(1.  -  . 005z ~ 1  -  .006z"2) (.000  -  . OOlz +  .500z~2) 
(1.  -  1 . 198z_1  +  .600zZ)(l  -  1 . 006z" 1  +  .401z"2) 


which  yields  the  numerator  polynomial 

z~2B(z)  =  .OOlz"3  . 500z * 4  -  .002z'5  -  .003z 


3.74 
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Further  analysis  shows  that  as  long  as  the  choice 
of  JL  is  no  more  than  two  units  less  than  the  number  of  units 


of  pure  delay  in  the  model,  the  cascade  structure  of  Figure  24 
allows  the  solution  z  H(z)  =  H(z)  regardless  of  system  order. 
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3.5  Parameters  of  TRAPS'  which  Provide  Enhanced 


Estimation  of  Seismic  Noise  Structure,  Using 

Arrays 

• 

The  frequency -wavenumber  (f-k)  power  spectral 
density  of  a  seismic  noise  field  is  one  common  way  to  summarize 

the  velocity  and  frequency  properties  of  the  noise  field  (39).  • 

The  determination  of  the  f-k  structure  of  seismic  waves  is  one 
of  the  most  powerful  techniques  available  for  the  study  of  the 
mode  structure  and  the  properties  of  these  waves.  Studies  of  • 

the  f-k  spectral  density  made  available  by  array  data  has  given 
insight  into  the  problem  (40-45).  This  section  focuses  on 

TRAPS'  Stage  One  noise  spectrum  and  the  f-k  delineation  made  • 

possible  by  time  tuning  specific  parameters  of  the  CAARMA 
model  of  the  preceeding  Sections. 

As  cited  in  Section  2.1  and  3.2  the  definition  of  • 

9 

noise  power  involves  a  number  of  assumptions  concerning  the 

time  and  space  stationarity  of  the  seismic  noise  field.  Also, 

each  event  is  different  to  some  extent  and  a  judgment  concern-  • 

ing  the  validity,  over  temporal  and  spatial  regions  of  interest, 

of  such  assumptions  must  be  made.  One  of  the  advantages  of 

TRAPS  is  that  due  to  the  adaptive  nature  of  the  processing  • 

system  most  assumptions  are  reduced  or  consequential  biases 

are  eliminated.  One  example  is  that  the  assumption  of 

stationarity  can  be  reduced  to  quasi-stationarity .  Also,  Stage  9 
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One  places  on  each  sensor  in  an  array  an  independent  adaptive 
filter.  These  filters  converge  to  those  data  values,  which 
were  made  local  to  that  sensor.  While  this  technique  has  far 
reaching  effects  in  yielding  an  increased  SNR,  this  independent 
adaptation  allows  automatic  removal  of  some  effects  of  seismo¬ 
meter  gain  variations.  For  the  filters  acting  on  the  individ¬ 
ual  geophones  of  the  array  a  sufficiently  long  adaptation  time 
was  selected  in  order  to  accomodate  the  variations  in  ambient 
noise  structure  without  disturbing  the  transient  signal 
signature.  For  the  filters  acting  on  the  beamformed  data, 
the  adaptation  time  varied  from  one  event  to  another  in  order 
to  incorporate  the  local  variance  of  each  event.  Furthermore, 
in  many  processing  techniques  a  too  small  of  a  window  or  model 
order  a  bias  is  injected  into  the  scheme.  If  the  length  is 
too  small  the  noise  will  appear  consistently  less  spatially 
organized  than  it  really  is.  While  the  order  of  the  CAARMA 
model  is  critical,  this  bias  is  bypassed  by  TRAPS.  This  is 
due  to  the  fact  that  the  adaptation  gives  the  CAARMA  model  the 
capability  of  a  model  with  much  greater  order  (46).  The 
adaptive  nature  of  TRAPS  also  eliminates  most  problems  asso¬ 
ciated  with  seasonal  changes. 


3.6  Technical  Results  of  CAARMA  model  on  NORESS  Data 


The  adaptive  ARMA  model  which  is  the  foundation  of 
the  CAARMA  model  and  the  slaved  first  Stage  which  prohibits 
the  effects  of  jamming  (nearby  man-made  disturbances,  local 
seismic  disturbances,  etc.)  and  signal  cancellation  each 
attacked  a  different  weakness  in  the  original  TRAPS.  Figures 
26  through  28  of  this  report  are  provided  to  allow  insight  into 
the  gain  which  each  of  these  two  improvements  provided.  The 
reader  is  referred  to  the  preceeding  semi-annual  report  if  it 
is  desired  to  obtain  a  separate  detailed  evaluation  of  each 
improvement,  since  only  the  results  of  the  complete  CAARMA 
model  shall  be  presented  at  this  time. 

The  complete  CAARMA  model  processing  system  enables 
TRAPS  to  detect,  locate  and  classify  seismic  events  for  nuclear 
monitoring  purposes  by  providing  the  necessary  information  for 
accuiate  results.  Since  the  Technical  Results  of  Section  4 
shall  be  demonstrated  by  means  of  TRAPS'  location  capability, 
this  Section  shall  present  the  results  of  this  Section  by  way 
of  instantaneous  power  spectral  plots.  Figure  29  is  a  raw 
data  plot  of  a  NORESS  event.  The  event  detail  is  listed  on 
the  Figure.  Further  event  information  can  also  be  found  in 
Table  2.  Figure  30. a  is  a  plot  of  the  instantaneous  power 
spectrum  of  the  event  represented  in  Figure  29,  using  the 
parameters  indicated  by  the  statistical  review  of  the  data  to 
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be  optimum  for  event  detection.  Figures  30.b-30.f  are  pro¬ 
vided  to  clarify  how  various  CAARMA  parameters  can  be  used 
to  pull  out  specific  features  of  an  event.  Since  the  results 
on  Section  4  shall  be  presented  in  a  manner  to  afford  proper 
evaluation  and  comparison  of  the  final  TRAPS,  Figure  30  is 


provided  so  that  visual  study  of  the  CAARMA  model  can  occur 
in  a  familiar  mode,  the  spectrum. 


4.0  TASK  3:  DETECTION,  LOCATION  AND  CHARACTERIZATION  OF 


NORESS  EVENTS  BY  TRAPS 


4.1  Introduction 


The  purpose  of  Task  3  for  the  contract  was  to  charac¬ 
terize  the  class  of  transients  which  compose  the  signal,  so 
that  automatic  identification  can  accurately  occur.  The  fre¬ 
quency-wave  vector  power  spectrum  describes  the  distribution 
of  power  versus  vector  and  temporal  frequency,  or  equivalently, 
versus  propagation  direction  and  temporal  frequency.  The 
frequency  wave  vector  spectrum  is  important,  since  the  detection 
and  bearing  estimation  of  the  source  is  based  upon  information 
contained  in  the  power  spectral  estimate.  Thus  among  the 
generic  signal  processing  concerns  of  an  array  system  is  the 
estimation  of  the  frequency-wavenumber  spectrum  of  the  signal 

9 

and  ambient  noise  field.  The  knowledge  of  these  spectrums  is 
essential  for  the  design  of  effective  detection  and  estimation 
algorithms,  when  the  received  signal  is  contaminated  by  ambient 
noise . 

In  order  to  demonstrate  the  logic  for  the  direction 
taken  by  the  Task  3  research,  the  equations  and  discussion  of 
Section  3  must  be  called  to  attention.  It  follows  from  these 
equations  that  the  CAARMA  coefficients  and  poles  and  zeros 
completely  determine  the  frequency-wavenumber  spectrum.  The 


0 


event  is  thus  represented,  not  by  a  set  of  discrete  values  of 
the  spectrum  or  its  integral,  but  by  a  number  of  parameters.  0 

These  parameters  being  the  poles  and  zeros  of  a  rational 
function  and  the  coefficients  of  the  CAARMA  model.  This  Sec¬ 
tion  demonstrates  how  this  information  is  utilized  to  yield  a  9 

highly  accurate  localization  algorithm  for  both  near  (events 
in  Sweden  and  Norway)  and  far  (Russian  events)  field  applica¬ 
tions.  Section  4.2  presents  a  method  for  improved  time  delay  * 

estimation  using  the  poles  of  a  CAARMA  model.  Section  4.3 
discusses  an  innovative  scheme  for  enhanced  detection,  location 
and  characterization  made  possible  by  the  coefficients  and  # 

residual  of  the  CAARMA  model.  Section  4.4  shows  the  results 
of  these  methods  on  NORESS/NORSAR  data. 

M 

$ 

* 

3 
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of  TRAPS  by  means  of  Pole/Zero  Decomposition 


The  semi-annual  report  described  a  localization 
algorithm  that  was  quite  successful  for  near  field  cases. 
Unfortunately,  the  time  delay  portion  of  the  algorithm  showed 
a  decline  in  performance  for  such  far  field  locations  as  Russia. 
The  semi-annual  report  demonstrated  how  the  poles  of  an  AR 
process,  due  to  their  signal  structure  sensitivity,  can  be  used 
to  detect  a  wave  front.  This  section  shall  expand  on  this  con¬ 
cept  to  illustrate  how  the  poles  and  zeros  of  an  AARMA  process 
can  be  used  to  accurately  determine  the  time  differences  of 
arrivals;  even  for  far  field  applications. 

A  variety  of  time  delay  estimation  techniques  have 
been  documented  which  adequately  address  the  problem  of  time 
delay  estimation  [47].  They  can  basically  be  classified  into 
two  categories.  One  is  the  crosscorrelation  method,  which  was 
the  method  originally  used  in  TRAPS  [48],  The  second  mode  of 
time  delay  models  the  time  delay  as  a  finite  impulse  response 
filter  with  one  process,  say  x(t),  as  the  input.  This  main¬ 
tains  similitude  to  an  MA  filter.  The  time  delay  estimate  is 
extracted  from  the  filter  coefficients.  The  innovational 
procedure  for  time  delay  estimation  to  now  be  described  is  one 
that  is  based  on  the  adaptive  ARMA  developed  under  this  contract 
and  detailed  in  earlier  sections. 


Before  the  correlation  is  established  between  the 
adaptive  ARMA  model  and  the  time  delay  analysis,  a  few  founda¬ 
tional  statements.  A  signal  is  assumed  at  all  geophones  within 
some  interval  of  time(T^,  .  (In  the  following  Section  it 
shall  be  demonstrated  that  this  assumption  is  not  required  for 
this  new  technique  to  work,  but  we  shall  maintain  the  assump¬ 
tion  at  this  time  in  order  to  facilitate  discussion.)  Let 
yQ(t)  be  the  filtered  time  series  measured  at  some  reference 
geophone.  Let  y^(t)  denote  the  filtered  series  on  the  second 
phone.  Then 


yo(t)  =  s  ( t )  +  nQ(t) 


yT(t)  =  s(t-Z  )  +  nx(t) 


where  t  (possibly  negative)  is  the  time  by  which  the  received 

signal  at  geophone  1  lags  the  received  signal  at  geophone  o  and 

s(t)  denotes  the  unknown  filtered  signal  waveform.  We  assume 

only  that  this  "signal"  is  within  the  interval.  nQlt)  and 

n^(t)  are  assumed  to  be  stationary,  uncorrelated,  zero  mean 

2  2 

Gaussian  processes  with  variances  ^  jnd  ,  respectively. 

\  nl 


The  time  argument  at  which  the  maximum  of  the 


function 
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h(t) 


&oyl  ,  ju)t  j 

sv7 Cw)e  dui 


occurs  is  equal  to  t.  This  estimator  is  equivalent  to  the 
Roth  processor  [49]  for  time  delay  estimation,  although  the 
attainment  is  different.  In  equation  (4.3)  Syoy^(c^)  is  the 
cross-power  spectrum  between  yQ(t)  and  y-^Ct)  and  Y0Y0  is  the 
power  spectra  of  yQ(t).  In  order  to  solve  t  in  ARMA  terms 


let  it  be  noted  that  [50] 


a,  e  -j“* 
k 


Sy  yi  v-  i  *  A 

5^r  CoJ)  =  £  "Sk'e  ’  b°  “  1  4*4 

k.-o 

Thus  the  time  delay  estimation  is  a  product  of  the  ARMA  process. 
This  method  was  shown  to  be  not  only  feasible,  but  relatively 
superior  to  most  methods,  as  long  as  the  SNR  was  greater  than 
OdB . 

Referencing  our  earlier  discussion  on  the  adaptive 
ARMA  model  of  TRAPS  and  the  section  that  follows,  the  time  delay 
problem  would  be  better  resolved  by  using  only  the  source  poles 
of  an  ARMA  model,  since  the  zeros  are  physical  related  to  the 
noise.  Also  means  must  be  available  for  multiple  source  time 
delay  estimation,  due  to  multipath.  In  1979  a  new  system 
identification  approach  for  time  delay  estimation  was  suggested 
by  Morf  [51],  using  a  multichannel  rational  model  for  the 
array  signals.  This  model,  which  is  useful  for  random  sources 
with  various  spectra,  was  later  developed  by  Porat  and 
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Friedlander  [52]  for  the  purpose  of  time  delay  estimation  of 
unknown  sources.  The  theory  was  based  on  an  ARMA  model.  Un¬ 
fortunately  not  only  did  amplitude  attenuation  need  to  be 
calculated,  but  the  receiver  noises  were  assumed  to  be  white 
and  spatially  uncorrelated. 

The  model  we  propose  to  use  in  TRAPS  for  time  delay 
estimation,  based  on  the  adaptive  ARMA  model  detailed  earlier, 
not  only  produces  the  time  delay  when  the  noise  may  be  corre¬ 
lated  among  the  sensors  and  in  time,  but  is  also  capable  of 
calculating  fractional  delays.  While  it  shall  be  demonstrated 
that  this  technique  is  not  only  highly  accurate,  it  must  be 
noted  that  an  advantage  of  the  technique  to  be  discussed  lies 
in  its  computational  simplicity.  Compared  to  the  generalized 
cross -correlation  technique,  the  new  time  delay  technique  does 
not  require  prior  knowledge  of  the  signal  and  noise  spectral 
densities.  The  new  method  also  has  the  advantage  that  it  is 
completely  independent  of  the  attenuation  factors. 

We  now  expand  on  equation  to  give  the  cross 

spectrum  (or  the  z-transform  of  the  cross  correlation)  between 
yQ  and  y^  in  a  multiple  source  environment  as  given  by 


Sy0yiU)  =  G  yQk  Gyjk  z 


‘  *yoylk 


(1/Z*) 


ak(z)ak*(l/z*) 


4.5 


where  ^yQy^k  =ZyQk  -  Zy^k  is  the  desired  time  delay  of 


arrival  of  a  signal  emitted  from  the  k  source  to  the 

receivers  y  and  y,  ,  and  Sn  ,  .. 

7o  7 1  *  J  yQy^  ^  represents  the  possible 

correlation  between  noise  components  T\yQ (t  Tg)  and  (t  Tg) . 

Here  t  denotes  the  time  index,  while  Tg  denotes  the  sampling 

interval.  Before  addressing  the  actual  time  delay  problem  it 

is  useful  to  expand  each  of  the  source  transfer  functions  into 

the  partial  fraction  decomposition 


.  >  C 


K  _  _  r-  _ 

av(z)  <L  1 

K  r-i 


^  kr 
^  kr  z  " ' 


where  \.Ckr}  are  the  residues  of  1  b^ (z) /a^ (z)}  at  the  corres¬ 
ponding  poles  [  A.  •  Substituting  4.6  into  4.5  now  gives 

Sy  yi(z)  -  IGyokGylk  z  Ay°yik^^C-~~ — z-i 

1  OK  IK  s._x  J  _^krZ 


C*ks  +  C  n 


t-Xk 


Multiplying  both  the  numerator  and  denominator  by  1  _A.kiz" 
and  considering  the  desired  limit  the  following  is  obtained 


A*  .^,y°ylU)  -  VZiyoylk 

3***1  S  ViO'3*) 


While  the  applications  of  these  equations  appear  conditional, 
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the  source  transfer  functions  may  have  repeated  poles  and/or 
be  nonproper.  This  can  be  proven  using  the  Laurent  expansion 
[53]  of  near  each  pole  by  L’Hopital's  rule.  The 

procedure  for  the  time  delay  estimation  results  from  equation 
4.8.  From  the  phase  part  of  equation  4.8  results  the 
following  time  delay  estimation 


a  y„y 


o7  lk 


»(-S  Wl  (1/1%!)}  -  0(5  wjjq)}  *  2tfr 


where 


)  -  tg'7-^r m) 


4.10 


The  ambiguity  term  in  4.9  is  resolved  by  the  amplitude  part  of 
4.8  which  is 


4Vik 


A  ^  yy 

log!6yoy1  u/rki>l  -  l0gts  yoy!  cjjeljj 

2  iog  |  ^  kl\ 
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4 . 3  Innovative  Acoustic  Signature  Identification  b) 

Means  of  CAARMA  Model  Output  for  Accurate  Location 
and  Classification 


The  purpose  of  this  Section  is  to  demonstrate  how 
the  information  made  available  by  the  CAARMA  model  of  TRAPS' 
Stage  One  and  Two  is  employed  to  yield  a  highly  accurate  local¬ 
ization  and  classification  algorithm  for  both  near  and  far 
field  applications.  The  procedure  begins  by  first  isolating  a 
seismic  event,  using  a  Weighted  Prediction  Residual  (WPR).  This 
WPR ,  which  will  be  presented  in  Section  4.3.1,  is  then  used  to 
elicit  the  signature  of  a  particular  event  with  the  aid  of  a 
segmentation  procedure.  This  WPR  segmentation  technique  is 
analyzed  in  Section  4.3.2.  The  WPR  views  the  waveform  from  a 
statistical  perspective,  while  the  segmentation  permits  the 
refinement  of  the  represented  waveform  a  structural  aspect. 

This  WPR  segmentation  procedure  is  able  to  classify  an  event  as 
an  earthquake  or  explosion,  use  large  phase  irregularities 
that  categorize  specific  field  data  to  improve  localization 
accuracy,  as  well  as  prepare  the  weighted  CAARMA  coefficients 
(WPC)  that  are  mathematically  presented  in  Section  4.3.1  and 
time  delay  information  provided  by  the  method  of  Section  4.2 
for  entry  into  a  feedback  system  with  a  maximum  likelihood 
scheme  for  an  extremely  efficacious  localization  algorithm 
that  will  be  detailed  in  Section  4.3.3. 
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4.3.1  Weighted  Prediction  Residuals  (WPR)  and  Coefficients 
(WPC)  Based  on  a  CAARMA  Model  to  Determine  Statis¬ 
tical  Waveform  Signature 

A  valid  supposition  is  that  such  statistical  ap¬ 
proaches  as  Residual  (LPR),  Linear  Prediction  Coefficients  (LPC) 
poles  and  zeros,  etc.  are  useful  tools  to  discriminate  and 
consequently  classify  an  event  (53-59).  A  superior  tool  in 
the  pattern  recognition  process  is  the  combination  of  a 
statistical  and  a  structural  approach  (60,  61).  As  a  result 
of  this  fact  information  provided  by  the  WPR  and  WPC  of  a 
CAARMA  model,  which  are  statistical  methods  based  on  the  LPR 
and  LPC  were  combined  with  a  structural  approach  to  accurately 
identify  the  several  waves  associated  with  an  event,  and  thus 
the  event.  So  that  the  combined  WPR/WPC  and  structural 
approach  may  best  be  understood,  the  logic  behind  a  WLPR  ap¬ 
proach,  resulting  from  a  CAARMA  model  as  well  as  the  mathema¬ 
tical  foundation  of  the  WPR  and  WPC  shall  be  established  in 
this  Section. 

The  basic  form  of  the  CAARMA  model  of  TRAPS  predicts 
the  windowed  signal  sample  x(n)  by  a  linear  combination 
of  $  excitation  and  (^response  time  series.  The  differential 
outcome  of  the  process,  which  in  turn  sets  the  foundation  for 
predicting  the  next  set  of  coefficients,  is  the  prediction 
error  e(n).  Recalling  the  equations  of  Section  3,  the 


adaptive  filter  constructs  a  sequence  of  prediction  coeffi¬ 
cients  a^  associated  with  the  excitation  and  associated 
with  the  response  time  series.  The  prediction  residual  is 
defined  as  the  sum  of  the  squares  of  the  sequence  e(n).  The 
magnitude  of  this  residual  is  both  a  function  of  the  prediction 
coefficients  {  a(n)  ,  b(n)^  ,  and  the  corresponding  excitation 
and  response  time  series.  For  a  given  data  sequence  Ix^n)} 
a  set  of  prediction  coefficients  {  a^  (n)V  can  be  calculated 
in  association  with  {b^(n)i  which  minimizes  the  residual, 
e  (n)  . 

A  method  quite  similar  to  this,  but  based  on  an 
AR  method  was  discussed  in  Section  5  of  the  semi-annual  report. 
It  was  quite  successful,  but  the  upgrade  to  the  CAARMA  model 
enhanced  the  method.  There  is  a  major  conceptual  difference 
between  the  prediction  residual  and  coefficients  of  an  AR  and 
a  CAARMA.  The  difference  lies  in  the  fact  that  the  residual 
and  coefficients  of  an  AR  model  involve  just  {a(k)\  , 
whereas  the  residual  and  coefficients  of  a  CAARMA  model 
involve  both  |_a(k)}  and  {b(k)}  .  The  problem  that  arises 
from  the  AR  to  CAARMA  transition  is  the  relevance  of  fb  (k)} 
in  the  task  of  waveform  classification.  The  relation  of  the 
CAARMA  coefficients  to  the  physical  laws  of  the  nuclear  moni¬ 
toring  problem  are  discussed  in  Section  3.2.  Once  compre¬ 
hension  of  the  logic  of  Section  3.2  occurs,  it  becomes 
evident  that  the  intent  of  ^b(k)]  is  to  reproduce  that 


portion  of  the  received  signal  which  is  the  noise.  Thus  it 
is  that  are  uniquely  related  to  an  event. 

Since  the  predictor  coefficients  fax(k)}  are 
uniquely  related  to  the  data  sequence  {x1  (n)}  ,  they  thus 
characterize  that  sequence.  That  is  the  coefficients  maintain 
parity  with  the  signature  of  the  event.  Now  suppose  another 
set  of  coefficients  {a 2  (k)}  are  calculated  by  the  CAARMA 
model  on  another  data  set  {x2(n)]  .  The  degree  of  similitude 

between  fa2Ck)l  and  (a^k))  determines  whether  or  not  the 
signatures  of  each  event  match.  Unfortunately,  much  of  the 


data  is  corrupted  by  artifacts  and  other  irregularities.  Thus 
despite  the  efficacy  of  the  CAARMA  model,  certain  coefficients 
and  residuals  can  be  corrupted.  The  resolution  to  this  prob¬ 
lem  is  the  appropriate  weighting  of  the  residuals  and  coeffi¬ 
cients.  Hence  the  terminology  WPR  and  WPC . 
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4.3.2  Segmentation  for  Structural  Def inition-Entropic 
Decision  Process 


The  method  of  Section  4.3.1  established  the  mode  by 
which  the  information  provided  by  the  CAARMA  model  can  be 
statistically  prepared  for  input  for  feed  back  into  the  range 
and  bearing  estimation  technique.  What  remains  to  be  disclosed 
is  the  mode  used  to  elicit  the  structural  information  of  the 
seismic  wave.  Seismic  waveforms  consist  in  segments  of  events 
of  unknown  duration,  starting  and  ending  points.  Hence  the 
technique,  now  presented,  is  based  on  the  segmentation  of  the 
WPR. 

The  time  trace  provided  by  an  array  of  geophones  is 
a  record  of  a  number  of  events  that  occur  after  source  energy 
is  sent  into  the  ground.  The  events  include  primary  reflec¬ 
tions  (compressional  and  shear,  multiples,  surface  waves,  etc.). 

♦ 

The  first  task  is  to  isolate  or  segment  these  so  that  dis¬ 
crimination  of  surface  waves  against  body  waves,  primaries 
against  multiples,  etc.  can  occur,  since  it  is  the  information 
revealed  by  each  of  these  events  that  is  necessary  to  yield 
accurate  range  and  bearing  estimates.  The  WPR  of  Section 
4.3.1  affords  the  segmentation  we  wish.  Figures  31  and  32 
give  a  visual  illustration  of  how  this  segmentation  occurs  in 
a  high  SNR  environment.  Performance  however  in  a  low  SNR  is 
sensitive  to  the  threshold  selected.  Some  trial  and  error 
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was  needed  to  obtain  a  proper  threshold.  Once  segmentation 
and  extraction  of  the  segment  transpired,  the  type  of  event 
now  needed  to  be  classified. 

So  that  each  type  of  event  is  accurately  classified, 
the  decision  process  which  minimizes  the  probability  of  in¬ 
correct  recognition  must  be  determined.  It  is  desirable  to 
begin  with  a  set  of  ideal  measurements  that  extract  only  those 
parameters  of  the  event  which  isolate  distinguishing  charac¬ 
teristics  of  the  event  type  identity.  Review  of  Section  3 
equations  and  the  discussion  of  Section  4.3.1  obviates  the  fact 
that  the  WPC  physically  represent  the  structure  of  the  earth 
and  consequently  the  distinguishing  characteristics  of  the 
event  type  identity.  A.  decision  was  thus  made  by  determining 
both  the  simple  average  of  the  WPR  of  an  event  with  a  reference 
WPR  and  the  entropic  distance  of  the  WPC. 

A  decision  by  a  simple  average  of  the  WPR  individ- 
ual  distances  made  a  comparison  between  the  reference  and  input 
contours  results  in  a  set  of  N  individual  distances  d^,  d^,..., 
d^,  then  the  overall  distance  D  is  given  by 


If  D  is  less  than  some  threshold  distance  Bt  then  the  decision 
is  made  to  accept  the  incoming  WPR  as  a  specific  type  of  event. 
Otherwise,  the  WPR  is  rejected.  Since  the  selection  of  contour 
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c 


measurements  and  corresponding  distances  is  based  to  a  large 

extent  on  intuition  the  selection  is  not  expected  to  be 

uniformly  optimal.  Therefore,  a  set  of  weights  w^,  w2,..., 

w  ,  were  used  such  that 
n 

N 

D  =  7*  w  d  4.13 

£=i  n  n 


where  D  best  discriminates  one  event  type  from  all  the  other 
event  types. 

As  for  the  entropic  distance,  let  a  denote  the  p- 
dimensional  column  vector  (a^ , . . . ,  a^)  and  b  denote  the  (p+1) 
dimensional  column  vector  (1,  a,,...,  a  ).  Now  consider  two 

x  P 

adjacent  sequences  of  measurement  for  reference  and  test  se¬ 
quences  XR  and  XT  of  length  NR  and  N^,  respectively.  The 
likelihood  function  under  the  null  hypothesis  that  the  two 
sequences  are  of  the  same  event  type  is 


Lo  =  ^p 


-frl 


R 


+  nt') 


exp[ -1/2 (Nr'  +  Nt')]  4.14 


P 


P 


where 


a  2 

Sp 


A  A 

=  a  1 C  a 
~p  -“P~P 
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and  the  poole  1  covariance  matrix 


is  given  by 
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4 


» 


0 


£p  =  (Nr'£r+Nt'£t)/(Nr'+Nt') .  4.16 

(62)  The  entropy  distance  d  can  then  be  defined  as 

L,  . 

d  -  -2  JU  (  °/L1)  4.17 


which  is  a  nonnegative  number  and  symmetrical.  If  cfR  ■  a"T 

A  \ 

and  olr=  ol^.,  then  d  =  0. 


>  .N 
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4.3.3  A  Maximum  Likelihood  Feedback  System 


As  cited  in  the  semi-annual  report,  the  original 
TRAPS  employed  a  maximum  likelihood  (ML)  method  to  solve  the 
location  problem.  Particulars  of  the  ML  mathematics  used  can 
be  found  in  Appendix  B.  The  method  of  maximum  likelihood  (63) 
was  introduced  by  the  English  mathematician,  R.A.  Fisher,  as 
a  general  approach  to  the  problem  of  estimating  the  "best" 
value  for  unknown  parameters,  0.  If  prior  knowledge  of  the 
probability  density  of  0  is  available,  the  likelihood  function 
should  be  multiplied  by  this  probability  density  before  the 
maximization  is  carried  out.  It  is  assumed  that  the  random 
parameter  is  restricted  to  a  range  of  values  and  that  the  pro¬ 
bability  density  of  the  parameter  is  constant  in  this  range. 

I'll 

The  arrival  time,  T^,  of  the  signal  at  the  i  geo¬ 
phone  can  be  expressed  by  the  following  equation: 

Ti  -  T  +  ~  [(xrx)2  +  (y-y)2]  l+  6i  4-l8 

where,  x^,  y^;  are  the  (x,y)  locations  of  the  i^  geophone 

x,y,T  are  the  location  and  origin  time  of  the  event 

V  is  the  effective  group  speed  of  the  signal 

is  an  additive  error  terms  assumed  to  be  zero 
mean,  Gaussian  with  known  covariance  matrix. 
Determination  of  the  arrival  times  at  each  geophone  by  the 


4-17 


method  of  Section  4.2,  allows  the  estimators  (x,y,T,V)  to  be 
chosen  so  as  to  maximize  the  likelihood. 


A  A 


L(x,t,T,V) 


(T  -  h.) 

A J  ^1 


A 


(T  -  hi) 


4.19 


where 


hi  5  T  +  v  t(xi  •  x)2  +  Cyi  ‘  y)25  1  4*20 

A  J 

However  if  A  =  (x,  t,  T,  v)  is  one  estimator  of 

A 

A 

this  sort,  and  B  is  another,  then  a  composite  minimum  variance 
unbiased  linear  estimate  C  can  be  obtained  as: 


A  a  . 

rv  A' 

J\.g  TV* 

A* 


4.21 


4.22  ’ 


Thereby  estimates  derived  from  the  independent  wave  components 
can  be  combined. 

The  investigation  of  this  technique  at  the  time  of 
the  semi  annual  report  revealed  this  method  was  to  have  a  gain 
significantly  in  excess  of  /”N  for  near  field  applications. 
Unfortunately,  for  far  field  cases  where  very  weak  events  were 
detected  much  of  the  signal  information  pertaining  to  propagation 
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effects  and  origin  was  lost  in  the  ambient  noise;  thus  pre¬ 
venting  Stage  Three  of  TRAPS  from  achieving  significant 
success.  What  was  apparent  was  that  the  TRAPS  process  must 
be  searched  to  understand  the  solution  to  this  dilemma.  Re¬ 
solution  came  from  two  basic  insights.  The  first  insight  was 
that  the  AARMA  model  contained  information  pertaining  to 
the  actual  physics  of  the  problem  and  effort  had  to  be  given 
to  improve  the  acquisition  of  the  information  to  insure  extreme 
efficacy.  Sections  prior  to  this  have  indicated  how  this  was 
accomplished.  The  second  insight  came  from  the  realization 
that  the  closed  loop  system  associated  with  the  generation  of 
the  adaptive  process  is  very  structured.  The  existence  of  this 
underly  ng  feedback  structure  was  extremely  useful  in  develop¬ 
ing  an  enhanced  TRAPS'  Stage  Three.  In  essence  what  these 
insights  allowed  us  to  do  was  to  first  extract  accurate  infor¬ 
mation  on  the  signal  even  in  a  low  SNR  environment  and  then 
use  a  feedback  system  to  isolate  the  parameters  that  allow  the 
ML  method  to  be  optimum,  e.g.  the  processing  system  was  in  an 
iterative  mode  with  the  physics  of  the  earth. 


4.4  Technical  Results 


The  results  of  this  Section  are  based  on  the  in¬ 
formation  made  available  by  the  CAARMA  model,  such  as  the 
WPR ,  WPC,  etc.  Thus  discussion  of  these  will  be  presented 
before  the  final  range  and  bearing  outcome  of  TRAPS. 

The  general  mode  of  the  WPR  method,  when  used  for 
event  signature  comparison  is  provided  in  Figures  31  and  32. 
Figure  32  is  actually  one  of  the  best  results  obtained  by  the 
LPR  method,  but  is  provided  since  it  graphically  demonstrates 
the  LPR  and  consequently  WPR  concept  so  extremely  well.  It  i 
important  to  recall  those  parameters  of  the  WPR,  which  have 
substantial  effect  on  the  results.  The  most  obvious  one  is 
that  of  the  hybrid  adaptive  model  order  and  adaptation  time 
used  to  generate  the  WPR.  Despite  the  spectrum  and  conse¬ 
quently  WPR  variability  generated  by  changing  the  order  and 
adaptation  time  of  the  adaptive  process,  the  variability  of 
the  WPP.  is  highly  consistent  for  each  order  group.  Thus,  as 
long  as  the  same  order  is  used  to  develop  all  the  LPR  classes 
for  pattern  recognition  comparison,  the  probability  of  iden¬ 
tification  remains  quite  accurate.  The  second  parameter  is 
that  of  the  number  of  residuals  of  which  the  WPR  is  summed. 
Whereas  this  factor  is  not  critical  in  that  there  is  a  sub¬ 
stantial  range  over  which  M  can  vary  before  results  deterio¬ 
rate,  it  is  still  an  important  issue.  In  general,  we  found 


the  summation  factor  of  three  times  the  order  of  the  CAARMA 
process  most  effective.  No  substantial  performance  deteriora¬ 
tion  occurred  between  1.0  and  50  times  the  order  of  the  CAARMA 
process . 

The  coefficients  of  the  CAARMA  model  are  affected 
by  the  same  factors  as  the  WPR.  Figures  33  and  34  are  plots 
of  the  data  that  was  processed  by  TRAPS  to  produce  the  co¬ 
efficients  found  in  Figures  35  and  36.  Figures  37  and  38  are 
these  same  coefficients  appropriately  weighted  for  anomoly 
removal.  These  WPC  and  the  associated  WPR  were  then  used  to 
perform  the  necessarv  range  and  bearing  calculations.  The 
range  and  bearing  results  can  be  found  as  follows: 

Table  3  -  Location  results  with  original  TRAPS  model 

4  -  Location  results  with  CAARMA  model 

5  -  Location  results  with  CAARMA  model  with 

improved  time  delay  of  Section  4.2 

6  -  Location  results  with  same  processing  scheme 

as  found  in  Table  5  with  Maximum  Likelihood 


feedback  scheme  added. 


5.0  TASK  4:  EXAMINE  TRAPS  FOR  INCORPORATION  INTO  A  COMPACT 


DIGITAL  SYSTEM 


When  the  signal  processing  research  efforts  progressed 
sufficiently  to  warrant  the  effort,  consideration  was  given  to 
design  features  for  software  and  hardware  requirements  of  a 
compact  digital  processing  system.  Since  the  realization  of 
such  a  system  within  a  real  time  limitation  was  possible  with 
the  original  TRAPS  (see  Figure  1)  (46)  ,  Task  4  efforts  thus 
needed  to  concentrate  solely  on  the  additional  software  and 
hardware  requirements  imposed  by  the  algorithm  improvements 
detailed  in  this  report.  Due  to  the  adaptive  foundation  of  the 
final  TRAPS  with  the  CAARMA  model  (see  Figure  3) ,  such  a  com¬ 
pact  system,  that  required  no  human  intervention  after  initial 
reference  waveform  setting,  was  totally  feasible.  Unfortunate¬ 
ly,  without  an  array  processor  as  part  of  the  system,  only 
near  real  time  performance  could  be  guaranteed. 


6 . 0  IMPLICATIONS  FOR  FUTURE  RESEARCH 


A  more  robust  TRAPS  I  for  automatic  detection,  location 
and  classification  is  possible.  This  could  be  accomplished 
first  by  means  of  research  including  Kalman  formalism  into 
TRAPS.  This  would  allow  the  continuous  upgrading  of  the  model, 
providing  time-varying  spectral  estimators  that  give  superior 
results  in  a  nonstationary  environment.  The  goal  of  a  more 
robust  TRAPS  I  could  also  be  met  by  incorporating  into  TRAPS 
a  Spatial  Dither  Algorithm,  which  adjusts  the  signal  to  remove 
distortion.  Since  the  signal  is  often  obscured  in  a  multi- 
path  environment  or  by  incomplete  data,  another  future  research 
project  could  address  these  problems  by  the  development  of  a 
theory  on  the  behavior  of  such  parameters  as  partial  correla¬ 
tion  coefficients,  pole/zero  migration,  etc.,  resulting  from 
TRAPS  I  and  the  relationship  of  the  behavior  to  the  signature 
of  the  source  of  interest  for  improvement  of  the  techniques 
presented  in  Section  4.3.  The  mathematical  theory  of  the 
supposition  proposed  is  based  on  the  coupling  of  the  theore¬ 
tical  concept  of  entropy  from  the  fields  of  Information  Theory 
and  Thermodynamics.  The  culmination  of  the  outcome  of  the 
above  stated  work  could  be  produced  onto  a  cognitive  imaging 
system,  such  that  TRAPS  could  be  given  the  ability  to  "educate" 
itself  on  required  input,  thereby  reducing  human  intervention. 
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Figure  4*  Instantaneous  Power  Spectrum,  resulting  from  TRAFS' 
Stage  Two  only,  using  a  filter  length  of  4*  This 
plot  was  the  product  of  a  simulation  study  specifi¬ 
cally  designed  for  a  THaHS'  Stage  One  processing 
gain  demonstration.  The  simulation  study  consisted 
of  3  seismic  waves  each  with  a  frequency  transition 
from  5.2  Hz  to  4*5  Hz  occurring  over  the  5  second 
length  of  the  signal.  Further  conditions  placed  on 
each  wave  were  as  follows: 

Signal  1:  5  db 

begins  at  second  17 
Signal  2:  10  db 

begins  at  second  h.2 
Signal  3:  5  db 

begins  at  second  71 

Superimposed  on  these  signals  were  three  .5  Hz,  10 
second  length  seismic  waves,  we  shall  term  noise, 
coming  from  the  same  direction  as  the  signals.  Fur¬ 
ther  conditions  placed  on  each  wave  were  as  follows: 
Noise  1:  15  db 

begins  at  second  17 
Noise  2:  40  db 

begins  at  second  42 
Noise  3:  10  db 

begins  at  second  71 


Figure  Instantaneous  Power  Spectrum,  resulting  from 
TRAPS'  Stage  One  and  Stage  Two,  using  a  fil¬ 
ter  length  of  4  for  each  stage.  This  plot 
was  the  product  of  a  simulation  study  speci¬ 
fically  designed  for  a  TRAPS'  Stage  One  pro¬ 
cessing  gain  demonstration.  The  simulation 
is  that  a3  detailed  on  Figure  4  caption. 
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6.  Condition  of  Power  Spectral  Density  (PSD)  Resolution 
for  Frequency  Resolution  Study;  the  results  of  which 

are  given  in  figures  7  and  27. 
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Figure  7.  Results  of  Adaptive  Beamformer  (AB)  Frequency  Resolution 
Evaluation,  Note  that  a  smaller  R  denotes  improved 
resolution. 
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Fourth  Order  Example:  Direct  Form  Algorithms. 
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Figure  19.  Fourth  Order  Example:  Cascade  Structure  Algorithms. 
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a-|.|  =  1 .5550  =  -.6490 

a12  ■  1.4990  a22  =  -.8480 


Adaptive  Parameters: 
a^  *  1 .4990  a21  » 
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Figure  16 


Fourth  Order  CRLMS  Example. 


Adaptive  Cascade  Structure. 


a )  Model : 

a1  ■  .0000  «2  ■  -.8000 

Adaptive  Parameters: 

•  -.0021  a2  ■  -.8010 

(After  300  samples) 

p  -  .01 


-2  0  2 

b)  Model : 

a1  -  1.2000  a2  ■  -.6000 

Adaptive  Parameters: 
a,  •  1.1997  a2  •  -.5998 

(After  800  samples) 

p  ■  .01 


-2  0  2 


c)  Model : 

»  -1.3000  a2  ■  -.4000 

Adaptive  Parameters: 
a,  -  1.2995  a2  •  -.3996 

(After  10,000  samples) 
p  »  .0005 


-2  0  2 


Figure  13 •  RLMS  Parameter  Tracks. 
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Figure  12.  RLMS  Parameter  Tracks. 
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Figure  11,  Performance  Surface  Examples. 
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Figure  2l\..  Modified  Adaptive  Cascade  Structure. 


# 


Data  Samples 
a)  l  »  4 


T 


T 


4A 


J _ l _ L_ 

5k  10k  15k 

Data  Samples 
b)  l  -  2 


Figure  2$.  Delay  Adjustment  Example. 


o  s  actual  zeroes 


/  «  actual  poles  s 
*  ■  calculated  zeroes 
+  »  calculated  poles 


*0/+  * 

I - 

0 


t  +/» 
H - 

5 

HZ 


e  +/ 

- - - Is 


Figure  26.  Adaptive  ARMA  pole-zero  estimates  compared  to 
actual  known  values,  after  adapting  to  3  data 
points. 
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Figure  27.  Results  of  Adaptive  Autoregressive  Moving  Average 
(ARMA)  Frequency  Resolution  Evaluation.  Note  that 
a  smaller  R  denotes  improved  resolution. 
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Figure  29a.  Display  of  Event  1  Received  Waveform 
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Figure  29c.  Display  of  Event  1  Received  Waveform 
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Figure  31,  Example  of  TRAPS-LPR  comparison  of  two  ambient  events 


Figure  32.  Example  of  TRAPS-LPR  comparison  of  two  similar  types  of  events 
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APPENDIX  A:  3-D  Hidden  Plot  program 
Apr  27  18153  1984  plt3d.f  Page  1 


c 

c... 

c 


10 

20 


c 


prog  ran  plt3d 
dimension  zmat  (91,  93) 
dimension  x(181),  z(lSl),  max 
data  param  /20.0,  30.0,  0,75> 
1  6.5,  6.75,  2.0/ 


\15>0;,  param  (V) 
0*0,  C.O,  0.0, 


initialize  variables 


deltax  =  param  (1) 
deltay  *  param  (2) 
deltaz  =  para*  (3) 
firstx  »  param  (4) 
firsty  ■  param  (5) 
firstz  *  param  (6) 
xlen  =  param  (7) 
ylmn  =  param  (8) 
zlen  =  param  (9) 
yincr  =  1.0 
res  =  40. 
n  =  91 

rotate  =  0.0 
angle  =  65.0 
iter  =  1 

do  10  i  =  1,  1500 

max  (i)  =  0 

continue 

do  20  i  =  1,  n 

x  < i >  3  float  (i  -  1) 

continue 

print  '('reading  data  please  wait')' 


c • • .open  input  file 
c 


open  <7,  file  =  'sl02',  status  -  > 

c 

c . . • read  all  data 
c 

\  read  (7,  1001) 

1001  format  (////). 
do  24  i  =  1,  91 

ict  s  o 

.jet  =  0 

do  22  j  »  1,  13 
ict  ■  jet  +  1 
.jet  «  ict  +  6 

read  (7,  1000)  (zmat  (i,  K),  K  =  ict,  .jet) 
22  continue 
24  continue 
1000  format  <7(f7.5,  lx)) 
c 

c...ring  bell 
c 

print 


c 

c»»»initiate  plot 
c 


9 


■^7 


r  ;-r  ■  w  ■-  ym 
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call  ginit 
call  clipon 

call  uvwlnd  (237.,  409?.,  o,, 
call  xywind  (0.0,  6.5,  0.0,  6.75.  \  r-) 
c 

c... clear  screen 
c 

call  erase 
c 

c .. .calculate  shifts  of  origin 
c 

rangle  =  angle.*  3.1415927  /  180,0 
yinc  =  yincr  /  deltay 
zshift  =  yinc 

xshift  =  yinc  *  cos  (rangle) 
c 

c » . .  rotate  plot 
c 

alph  =  rotate  *  3.141527  /  ISO  0 
sina  =  sin  (alph) 
cosa  =  cos  (alph> 
c 

c.t.main  loop  to  plot 
c 

do  120  K  -  1,  n 
c 

c..,take  data  one  line  at  a  time 
c 

iflag  «  1 
do  110  .j  =  1,  n 
z(j)  =  zmat  (.j,  h) 

>!*  “  “  first::)  /  delta  ;  •)  -ioat  ( iter  -  1)  %  >; shift 

^1  -  (z:.j)  -  firstz)  f  deltaz  +  fie  at  (;te-'  -  '  )  *  z'r,‘ii~t 

if  (.j  .no.  1)  qr,  to  30 
x2  =  xl 
z2  =  zl 

>:2p  *  cosa  *  >:2  -  sina  *  z2 

z2p  =  cosa  *  z2  +  sina  *  x2 

>:9  *  x2p 
z9  =  z2p 

call  xymove  (x9,  z9) 
c 

c ...  interpolate,  pick  integral  endpoint* 
c 

slope  =  0.0 
go  to  40 
30  continue 

slope  =  (zl-z2)  /  <xl-:;2) 

40  continue 

b  *  zl  -  (slope  *  >; l .) 
nstart  =  int  (x2  *  res) 
nend  -  int  (>:1  *  res) 

>:2p  =  cosa  *  x2  -  sina  *  z2 
z2p  =  cosa  *  z2  +  sina  *  >;2 
c 

c...plot  at  points  corresponding  to  values  in  max  army, 


# 


i 


l 


■■■  , 

«v»s  iM 
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c.«,as,¥#ll  a*  x2 1  if  visible 

■  •  f*'  K  <  ;*  fy 

'■yi'-  toyi oo  i  *»  nstart,  nend 

*tm  float  <i)  /  res 
*3  ■  slope  »  x3  +  b 
*3p  *  cosa  t  x3  -  sina 
t3p  *  cosa  *  z3  +  sina 
iz4  ■  int  (z3  *  400.0) 
if  (nstart  ,ne.  i  .or, 
i*2hot  ■  max  (i)  -  int 
iz5  ■  int  (z2  *  400.0) 


*  -  v  -  *.  *  i.*  ••  .;*•  . 


23 

x3 


i  .eq.  1)  go  to  60 
(res  $  (max  <i)  -  aiax 


, 

if  {iz5  .It.  iz2hat)  go 

caliinydro*  (x2p, 

z2p) 

r 

iflCf  «  0 
go  to  60 

50 

continue 

.  A".- '  - 

call  xymove  <x2p, 

z2p> 

iflag  =  1 

60 

continue 

if  (iz4  .It.  max  <i))  go  to  30 
if  (iflag  .eq .  1)  go  to  70 
coll  xydraw  <x3p,  z3p) 

•ax  (i)  =  iz4 
go  to  100 
70  continue 
iflag  =  0 
go  to  90 
B0  continue 
iflag  =  l 
90  continue 

call  xymove  <x3p,  z3p ) 

100  continue 


c  .  •  .update  x2,  z2»  iteration  counter 
c 

22  =  2 1 
x2  =  x.1 

110  continue 

iter  =  iter  +  1 
120  continue 


c 

c... alpha  mode  and  ring  bell 
c 

call  xymove  (0.0,  6.75) 
call  size  (2) 
call  alpha 
print  ' ( ** ) ' 
c 

stop 

end 


(i-1))  *  <x3-x2) ) 
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Stage  Three  of  TRAPS  Filter:  Range  and  Bearing 

Estimation 

The  adaptive  filter  accurately  determined  the  existence 
of  the  signal  within  some  time  interval.  This  time  interval  was 
then  used  as  an  input  for  time  delay  estimation.  The  exactness 
of  the  adaptive  filter  in  specifying  the  signal  permitted  the 
time  delay  procedure  to  reliably  specify  the  relative  arrival 
times  at  each  acoustic  sensor.  Once  the  arrival  times  of  the 
signal  at  each  sensor  were  estimated,  a  maximum  likelihood 
method  was  able  to  estimate  the  range  and  bearing  of  the  source. 
We  will  detail  the  mathematics  at  this  time. 

As  briefed  in  the  above  paragraph,  the  first  two  stages  of 
*  TRAPS  allows  one  to  stipulate  the  condition  of  a  signal  at  all 
geophones  within  some  interval  of  time  (Tj_,  T2)  .  We  can  then 
estimate  their  relative  arrival  times  by  pairwise  cross  correla¬ 
tion  of  a  filtered  series  measured  at  each  geophone  with  that 
of  a  common,  arbitrarily  chosen,  reference  phone.  Let  Y0(t) 
be  the  filtered  time  series  measured  at  the  reference  geophone, 
and  let  Y-^t)  denote  the  filtered  series  on  the  second  phone. 
Then 

Y0(t)  =  S  (t)  +  nQ(t) 

and 

Y^t)  =  S(t-x)  +  nx(t) 

where  t  (possibly  negative)  is  the  time  by  which  the  received 
signal  at  geophone  1  lags  that  at  geophone  0,  and  S(t)  denotes 
the  unknown  filtered  signal  waveform.  We  assume  only  that  this 
"signal”  lies  within  the  interval.  More  precisely,  the  filtered 
wave  form  is  assumed  to  be  similar  on  each  geophone;  some  dissi¬ 
milar  residue  in  the  signal  may  be  combined  with  the  noise  so 


A-l 


long  as  it  does  not  greatly  affect  our  assumptions  regarding 
the  noise.  We  assume  the  noise  nQ(t)  and  n^(t)  to  be  stationary, 
to  be  uncorrelated,  and  to  be  zero  mean  Gaussian  processes  with 
variances  o^o  and  respectively. 

We  choose  the  delay  x  which  maximizes  the  a  postori 
probability  density  of  t  given  the  measured  series  yQ,  y^. 

That  is,  we  choose  x  such  that  p(x/yQ,  y^)  is  maximized. 

It  can  be  shown  (14)  by  this  criterion  that  A  »  t  -  t*, 
the  estimated  minus  the  true  delay,  has  the  standard  deviation 


where  W  is  the  rms  bandwidth  in  Hz  and  — - —  is  twice  the 

N 

signal  energy  to  noise  energy  ratio.  For  the  NORESS  data 
listed  in  Table  1,  we  achieved  an  average  standard  deviation 
of 


a.  -  .018  sec. 
aAVE 

This  cross  correlation  procedure  requires  that  the  signal 
be  known  to  lie  within  some  interval.  The  adaptive  filter 
detailed  in  Sections  2.2  and  2.3  was  able  to  accurately 
identify  the  interval  in  which  the  signal  existed. 

A  set  of  relative  arrival  times  has  now  been  obtained. 
These  times  are  now  processed  by  the  third  stage  to  determine 
if  they  jointly  indicate  that  their  respective  transients 
all  came  from  the  same  location  in  space.  The  criterion  as 
to  whether  this  is,  in  fact,  the  case  is  given  by  the  area 
of  the  x,y  confidence  ellipse.  The  definition  of  the  x»y 
confidence  ellipse  is  that  elliptic  contour  within  which  we 
can  be  95%  certain  that  the  transient  wave  front  originated. 


The  equations  that  follow  explicitly  state  how  the 
x,y  position  is  obtained  from  the  relative  arrival  times, 
and  how  the  confidence  ellipse  is  determined.  Also 
presented  is  the  statistical  means  by  which  the  various 
geophone  channels  are  combined  to  obtain  estimates  which 
are  more  precise.  The  times  T^  are  those  obtained  by 
the  cross  correlation  procedure  just  described.  Note 
that  the  time  T^  to  the  reference  phone  is  zero. 

The  fundamental  physical  model  for  the  original  TRAPS  is 
both  simple  and  effective.  The  program  is  not  however,  in  any 
way  commited  to  these  equations  and,  as  will  be  discussed  in  later 
sections,  these  equations  were  made  more  sophisticated  by  the  re¬ 
search  efforts  of  this  contract. 

As  depicted  in  Figure  suppose  a  source  originating 

at  time  T  is  located  at  position  x,y  and  its  seismic  wave 
is  received  at  times  T^ :  i  =  1,...,N.  Let  the  ranges  to 
these  sensors  be  denoted  by  r^ :  i  =  1,...,N.  Then  to 
first  order  (assuming  straight  line  travel) ,  we  have  the 
following  equations: 


Ti  =  T  +  ^  jjxj-x)  2  +  (yi-y)2 

These  comprise  n  equations  in  4  unknowns  (x,  y,  t,  v) ,  and 
if  N=4  they  can  be  solved  exactly.  However,  the  arrival 
times  in  fact  contain  random  noise  perturbations  which  will 
therefore  introduce  error  into  the  solution.  Consequently, 
many  observations  are  made  and  the  solution  which  maximizes 
the  likelihood  function  is  chosen. 

The  model  this  becomes: 


T^  =  T  + 


(Xi-x)2  +  (yi-y)2 1  +  ti 


where  the  Gaussian  random  process  e  has  zero  mean  and  known 
covariance  matrix  A.- 


0 


We  choose  x,  y,  T,  v  to  minimize  the  form 


Q(x,y,T,v 


>  -  tf,  (Ti-hi>[Aj1]  u  ‘W 


where 


“i  "  T  +  V  [ 

A; 


(x^x)2  +  (yi~y)2J 


i  = 


denotes  the  inverse  of  .  ,  _ 
6  J  £ 


A, 


We  now  let:  x  *  (x,y,T,v) 


3hi 


3xj 


The  normal  equations  for  this  problem  become 


N 

‘if 


Ar 


ij 


(Tj-hj)  =  0 


k  =  1, .  . .  ,  4 


These  non-linear  equations  are  to  be  solved  for  the  unknown 
vector  x  estimating  x. 

The  solution  is  obtained  as  the  convergent  limit  of  the 
Gauss-Newton  iteration  scheme 


x  (n+1)  =  ;<n)  + 


[(hTA  A  H)'lHTA;1]Mn)£(n) 


where  x^  denotes  the  value  of  x  at  the  nfc^  iteration 


r(n)  -  T  -  h  (x(n)) 


»-4 


9 


T  -  (T^#T2# • • • *Tn) 


£  is  an  initial  guess 


The  iterative  scheme  will  converge  for  virtually  any 
choice  of  to  the  same  limit.  In  the  program  implemented 

the  matrix  is  taken  to  be  oj^g  times  the  identity  matrix. 


-  |  T±  -  hi  (x(n))  2 


The  estimator  jc  has  a  covariance  matrix  given  by 


Aa 


<«tA; 


The  lo  confidence  ellipse  about  x  is  derived  from  this  matrix. 

The  above  method  estimates  a  velocity  which  should  be 
taken  to  represent  the  mean  horizontal  group  velocity  from 
the  source  to  the  array.  To  a  great  extent  the  method 
succeeds  on  real  data  because  this  mean  velocity  is  a 
very  stable  physical  parameter,  and  its  value  is  not  assumed 
to  be  known  bv  the  p roar am. 

The  equations  governing  the  combination  of  observations 
are  as  follows: 

Let  x  and  y  independently  estimate  a  parameter  a.  Let 

and  *  denote  their  covariance  matrices.  Then  the  minimum 

y  a 

unbiased  composite  linear  estimator  Z  is  given  by 


Az  •  i  -  A,  •  » +  A;  •  x 


APPENDIX  C 

INFORMATION  RELATED  TO 

NORESS  SEISMIC  ARRAY  DATA 

.  '  • 
t  * 

Seismic  data  has  been  collected  in  one  of  two  array 
configurations: 

a.  A  5  element,  3  axis  per  element  array  during  the 
time  period  10  June  1983  to  5  July  1983. 

b.  A  21  element,  vertical  axis  array  starting 
approximately  1  September  1983. 

The  data  sample  rate  is  40  Hz  and  the  system  displacement 
response  is  depicted  by  Figure  1. 

The  three  axis  systems  (Z,  N,  E)  were  located  at 
positions  1,  2,  6,  9  and  11  used  in  original  NORESS 
configuration  (see  Figure  2) .  Locations  are  given 
referenced  to  position  #2.  Position  #2  is  located  at 
60.73523°  North  latitude  and  11.54136°  East 
longitude.  Location  data  are  provided  in  Table  I.  In 
Table  I  ana  Table  II,  a  positive  X  coordinate  is  East  of 
the  reference  point  and  a  positive  Y  coordinate  is  north 
of  the  reference  point. 


# 


-2- 

e 

The  21  element  vertical  axis  array  was  deployed  at  points 
near  the  points  to  be  used  in  the  planned  DARPA/DOE/ 
NORSA&  1984  array.  (See  Figure  3).  The  center  location 
(A.  or  Channel  1)  is  near  the  reference  location 
(position  #2)  defined  for  the  3  axis  array.  Element 
locations,  referenced  to  the  center  location  are  defined 
in  Table  11. 


» 
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TABLE  I 

3  Axis  Array  Locations 


Channel 

3  3 

NORESS 

Location 

T>A*IA  “TAPf 
CUAHKKlL  xd 

Axis  - -  X  (EW) 

Y  (NS 

1 

2 

Z 

I 

A2  0 

0 

2 

2 

N 

ZN 

3 

2 

E 

ZE 

4 

1 

Z 

IX  -62 

+271 

5 

1 

N 

1  A/ 

6 

1 

E 

IS 

7 

8 
5 


6 

6 

6 
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TABLE  II 

VerCical  Axis  Array  LocaCions 

P47A  T Aft  » 

CHAHW-  XV> 

Channel  NORESS  Axis  — Location  (Meters) 

# _ Location _ X(EW) _ Y(NS) _ Z 


0 


1 

Center  (A  >  Z 

o 

A* 

0 

0 

0 

2 

A1 

Al 

18.2 

144.6 

-15.3 

3 

a2 

*1 

124.4 

-57.1 

1.3 

4 

*3 

-144.9 

-66.4 

“8.1  -/ 

5 

*1 

81 

63.5 

270.5 

-2.7 

6 

“2 

81 

316.0 

-7.6 

11.5 

7 

-3 

83 

204.9 

-234.4 

15.5 

8 

84 

-185.3 

-245.2 

-1  -V 

9 

*5 

8r 

-266.1 

127.0 

-12.3  : 

10 

C1 

C.I 

59.0 

693.5 

9.7  :.;5 

11 

C2 

C2 

623.8 

361.0 

39.7  :  ■ :- 

12 

C3 

C3 

658.3 

-207.9 

56.5 ^ 

13 

=4 

c4 

202.1 

-671.6 

6.8  . 

14 

"3 

cr 

-360.6 

-624.6 

-4.1 

15 

C6 

-770.5 

-82.9 

io.o 

> 

16 

C7 

Cl 

-454.9 

571.2 

-26.4  9  :~~A 

17 

w 

P 1 

138.4 

1476.4 

5.9 

18 

“2 

01 

1060.2 

1100.5 

67.6 

19 

“4 

P4 

1225.1 

-933.4 

81.6 

20 

»s 

or 

-412.2 

-1466.6 

45.3  , 

21 

d7  z 

P7 

-1347.2 

-734.7 

32.3 

22 

(Wind  Speed) 

vUS 

0 


0 


0 
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Table  1.  Periodogram  Resolution 
as  a  Function  of  SNR 


MEAN  OF 

STANDARD  DEVIATION 

*n 

0F  *N 

00  0.031374  ±0.0 

20  dB  0.031374  ±0.000082 

10  dB  0.031359  ±0.000258 

0  dB  0.031328  ±0.000798 


-10  dB 


0.030989 


±0.002362 


Table  3.  Original  TRaPS  Bearing  Results 


Event  No. 

1 

2 1 

3 1 

V  - 

.5* 

6* 

7* 

3* 

Azimuth 

Accuracy 

tl° 

r.05° 

il° 

11° 

ii° 

±3° 

1 3° 

o 

U\ 

O 

• 

♦i 

Table  1^..  TRAPS  Bearing  Results  with  CAARMa  Model 


Event  No. 

1 

2' 

3' 

it* 

5* 

6« 

7* 

8» 

Azimuth 

Accuracy 

1 1° 

±.05° 

tl° 

i  1° 

1 1° 

+  1° 

tl° 

i.05° 

Table  5.  TRaPS  Bearing  Results  with  CaaRMa  Model 
and  Improved  Time  Delay  Method 


Event  No. 

1 

2' 

3' 

6' 

7' 

8* 

Azimuth 

Accuracy 

±  .05c 

o 

U\ 

O 

. 

t  1° 

tl° 

tl° 

i  1° 

i  1° 

Hi 

Table  6.  TRAPS  Bearing  Results  with  CAaRMa  Model, 
Improved  Time  Delay  &  ML  Feedback  Scheme 


i 

Event  No. 

1 

2« 

3' 

k' 

5« 

6» 

7'  _ 

3 » 

azimuth 

Accuracy 

RS 

Hi 

IPS 

HB 

ms 

bp 

m 

